93-00102 


AFIT  /GEP/ENP/92D- 1 


DTIC 


ELECTE 
JAN  7  1993 


A  TWO-DIMENSIONAL  PARTICLE 
SIMULATION  OF  PARALLEL  PLATE 
RADIO-FREQUENCY  (RF) 
GLOW  DISCHARGES 

THESIS 

Eric  J.  Bennett,  Captain,  USAF 
AFIT/GEP/ENP/92D- 1 


Approved  for  public  release;  distribution  unlimited 


93  1  04  I. 


AFIT/GEP/ENP/92D- 1 


A  TWO-DIMENSIONAL  PARTICLE  SIMULATION  OF  PARALLEL  PLATE 
RADIO-FREQUENCY  (RF)  GLOW  DISCHARGES 


THESIS 


Presented  to  the  Faculty  of  the  School  of  Engineering 
of  the  Air  Force  Institute  of  Technology 
Air  University 
In  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 
Master  of  Science  in  Physics 


D^QD/.LIr:-  .... 


Eric  J.  Bennett,  B.A.,  M.S. 

Captain,  USAF 

December  1992 

Approved  for  public  release;  distribution  unlimited 


Ac«««3l(M  For 

NT  IS 

*?4<5  fli> 

UaAnuoun«u4 
J ua  t  i  i'  *.* 


n 

□ 


By -  - 

Di  •  !  s.j  r  i  $:/ 
Av.»i  .  i;  •  r  t- 


;di  st 

k-\ 


.  1  «s 


•  J  r 


'I 


1 


Preface 


The  purpose  of  this  research  was  to  develop  a  two-dimensional  particle  simulation  of 
the  parallel  plate  RF  glow  discharges  used  by  industry  in  plasma  assisted  processing.  The 
need  for  this  work  arises  from  the  recognized  lack  of  understanding  regarding  the 
relationships  between  the  external  operating  parameters  and  internal  processes  in  the  RF 
glow  discharge  (Shohet,  1991:726).  The  National  Research  Council  has  recommended 
that  a  national  initiative  be  undertaken  to  increase  the  level  of  scientific  knowledge 
regarding  RF  glow  discharge  reactors.  In  particular,  the  National  Research  Council  has 
recommended  that  multi-dimensional  models  be  developed  to  study  the  RF  glow 
discharge  and  aid  in  future  reactor  designs  (National  Research  Council,  1991:51). 

Since  the  calculation  is  computationally  intensive,  the  initial  model  was  developed 
with  the  mind  set  of  including  only  those  processes  essential  for  preliminary  testing. 

Initial  testing  indicates  that  the  model  reproduces  many  of  the  qualitative  features 
predicted  by  some  of  the  existing  one-dimensional  models  and  theory.  Further  work  is 
required  to  gain  a  more  quantitative  agreement.  As  it  now  stands,  the  model  runs  on  a 
computer  which  limits  the  duration  of  a  simulation  and  the  number  of  particles  in  the 
simulation  (due  to  memory  constraints).  Ultimately,  the  model  should  be  optimized  for 
run-time  and  ported  to  a  super  computer  class  machine  for  computer  experiments 
consisting  of  longer  run  times  and  a  larger  number  of  particles. 

Chapter  one  provides  background  on  the  RF  glow  discharge  and  motivates  this  work 
by  outlining  industrial  dependence  on  plasma  assisted  processing.  A  specific  example  in 
the  microelectronics  industry  is  used.  Chapter  two  outlines  the  various  general  methods 
available  for  RF  glow  discharge  simulation  and  reviews  some  previous  models.  In 
addition,  the  simulation  method  used  in  this  research  (particle  simulation)  is  discussed  in 
detail  and  the  specific  algorithms  used  in  the  model  are  described.  Chapter  three  provides 


u 


the  results  of  some  simple  test  cases  which  give  confidence  in  the  proper  operation  of 
portions  of  the  model.  Also,  the  results  of  initial  experiments  conducted  using  the 
simulation  are  discussed.  Chapter  four  provides  required  changes  to  the  model  and 
suggestions  for  future  work. 

I  would  like  to  thank  my  thesis  advisor  Dr.  W.  F.  Bailey  for  introducing  me  to  this 
research  topic  and  for  his  encouragement  and  assistance.  I  would  also  like  to  thank  my 
committee  members.  Dr  A.  Garscadden  and  Capt  M.  Stoecker,  for  their  suggestions  and 
comments  on  this  document.  As  always,  much  gratitude  goes  to  my  wife  Cindy  for  her 
patience  and  understanding. 

Eric  J.  Bennett 
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Abstract 

A  two  dimensional  model  of  parallel  plate  RF  glow  discharges  was  developed  to 
study  discharge  phenomena  important  in  plasma  assisted  processing  of  materials.  The 
particle-in-cell  method  is  used  to  calculate  the  trajectories  of  computer  particles  under  the 
influence  of  both  self  and  applied  fields.  Monte  Carlo  methods  using  the  null  collision 
technique  are  used  to  model  collisions  between  charged  particles  and  neutral  gas  atoms. 
Results  of  computer  experiments  are  presented  with  special  emphasis  placed  on  ion 
motion  in  the  sheath  regions.  Experimental  results  show  some  qualitative  agreement  with 
one-dimensional  model  results.  Further  work  required  to  gain  quantitative  agreement  is 
outlined. 


A  TWO-DIMENSIONAL  PARTICLE  SIMULATION  OF  PARALLEL  PLATE 
RADIO-FREQUENCY  (RF)  GLOW  DISCHARGES 


I.  Background  and  Motivation 

Beyond  pure  scientific  curiosity,  why  study  the  RF  glow  discharge?  The  answer  is  - 
major  manufacturing  industries  have  become  increasingly  dependent  upon  the  RF  glow 
discharge  for  a  wide  variety  of  materials  processing  techniques  (National  Research 
Council,  1991:4-5).  Furthermore,  basic  scientific  knowledge  of  the  discharge  is  seriously 
lagging  the  need  for  technology  development.  This  chapter  introduces  the  physical 
properties  of  the  RF  glow  discharge  and  establishes  the  growing  dependence  of  industry 
on  plasma  processing,  as  well  as  the  need  for  increased  scientific  understanding  of  RF 
glow  discharge  phenomena  (specifically,  the  requirement  for  multi-dimensional  computer 
models). 

Characteristics  of  the  RF  Glow  Discharge 

The  typical  parallel  plate  RF  glow  discharge  is  shown  in  figure  1.1  (Flamm  and  Herb, 
1989: 14-17).  Flamm  and  Herb  provide  a  good  overview  of  some  of  the  major 
characteristics  of  the  RF  glow  discharge.  A  gas  at  low  pressure  (typically  ten  mtorr  to 
one  torr)  is  introduced  into  the  space  between  two  parallel  plate  electrodes.  A  plasma  is 
produced  via  the  application  of  RF  power  through  the  electrodes.  The  electrodes  are 
typically  driven  at  the  industry  standard  frequency  of  13.56  MHz  and  a  peak  voltage  of  a 
few  tens  to  one  thousand  volts.  The  plasma  consists  of  a  weakly  ionized  gas  with 
approximately  one  charged  particle  for  every  105  -  106  neutrals.  The  applied  RF  power  is 
the  source  of  energy  for  the  production  of  charged  particles.  The  charged  species  consist 
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of  heavy,  relatively  immobile  positive  ions  (the  positive  charge  carriers)  and  light,  highly 
mobile  electrons  (the  negative  charge  carriers),  although  some  plasmas  contain  negative 
ions. 

Flamm  and  Herb  also  describe  the  two  distinct  spatial  regions  which  form  in  the  RF 
glow  discharge  (Flamm  and  Herb,  1989: 15-17).  The  majority  of  the  discharge  is 
electrically  neutral  (the  plasma  bulk),  consisting  of  approximately  equal  numbers  of 
positive  and  negative  charges  with  a  number  density  of  10 15  - 1018  nr3.  The  other  region 
of  the  discharge  is  formed  between  the  bulk  plasma  and  any  surface  which  comes  into 
contact  with  the  plasma  such  as  the  walls  of  the  discharge  tube  and  the  electrodes.  These 
are  known  as  sheath  regions.  The  formation  of  the  sheath  which  forms  between  the 
plasma  bulk  and  the  radially  bounding  glass  wall  can  be  described  as  follows:  in  the 
sheath  region  diffusion  of  charged  particles  to  the  bounding  surface  severely  reduces  the 
charged  particle  densities.  Since  the  electrons  diffuse  faster  than  the  ions  (owing  to  their 
much  smaller  mass  and  higher  mobility),  a  net  positive  space  charge  region  is  left  behind 
which  forms  an  electric  field  that  promotes  ion  diffusion  and  retards  further  electron 
diffusion.  Eventually  the  electrons  and  ions  diffuse  to  the  wall  at  an  equal  rate.  The 
potential  difference  associated  with  the  radial  sheath  is  much  smaller  than  that  of  the  axial 
sheath  which  forms  between  the  electrodes  and  the  plasma  bulk.  Due  to  the  depletion  of 
electrons  in  the  sheath  regions,  the  sheaths  are  poor  conductors  as  compared  to  the  charge 
rich  plasma  bulk.  This  situation  leads  to  large  potential  differences  across  the  axial 
sheaths.  In  fact,  "most  of  the  potential  drop  appears  across  the  sheath"  (Flamm  and 
Herb,  1989: 17).  Since  the  sheath  region  is  narrow,  it  is  characterized  by  high  electric 
field  strengths  as  compared  to  the  bulk. 

The  axial  sheath  regions  play  an  important  role  in  the  RF  glow  discharge.  The 
sheaths  that  form  between  the  electrodes  and  the  plasma  bulk  in  an  RF  glow  discharge 
vanish  and  reform  once  during  each  RF  cycle  due  to  the  application  of  a  time  varying 
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voltage  to  one  of  the  electrodes.  This  sheath  motion  imparts  energy  to  the  charged 
particles.  Since  the  ions  are  much  heavier  than  the  electrons,  they  are  unable  to  respond 
to  the  instantaneous  potential  at  high  driving  frequencies  —  the  heavy  ions  respond  to  a 
time-averaged  potential.  Thermal  ions  entering  the  sheath  are  accelerated  by  the  time- 
average  sheath  electric  field  towards  the  electrode  until  they  consequently  strike  the 
electrode.  The  lighter  electrons  respond  much  differently  to  the  time-varying  sheath 
which  imparts  large  amounts  of  energy  to  the  electrons. 

Boeuf  and  Belenguer  explain  that  electrons  can  gain  energy  from  the  plasma  sheath 
by  any  of  three  mechanisms  (Boeuf  and  Belenguer,  1990: 164  -  166).  Depending  upon 
the  discharge  conditions  any  one  of  these  mechanisms  may  be  the  dominant  one.  One 
way  in  which  electrons  gain  energy  from  the  sheath  is  by  secondary  electron  emission 
due  to  ion  bombardment  (Boeuf  and  Belenguer,  1990:  166).  When  an  ion  strikes  an 
electrode  a  secondary  electron  may  be  emitted.  This  secondary  electron  created  at  the 
surface  of  the  electrode  is  accelerated  through  the  large  sheath  field  away  from  the 
electrode  and  towards  the  central  region  of  the  plasma.  These  secondaries  enter  the  bulk 
with  high  energies. 

A  second  manner  in  which  electrons  gain  energy  is  through  what  Kushner  terms 
"wave  riding"  (Kushner,  1986:  188).  As  described  above,  during  a  portion  of  the  RF 
period  the  sheath  contracts,  allowing  some  higher  energy  bulk  electrons  to  gain  access  to 
the  sheath  region  where  they  are  subsequently  pushed  back  into  the  bulk  as  the  sheath 
reforms.  The  electrons  that  intrude  upon  the  sheath  region  may  reach  high  energies  as 
they  "surf"  on  the  expanding  sheath  back  to  the  bulk  (Boeuf  and  Belenguer,  1990:  166). 

Lastly,  the  sheath  imparts  energy  to  the  electrons  through  electron-sheath  collisions 
(Boeuf  and  Belenguer,  1990:  166).  As  the  sheath  expands  away  from  the  electrode  the 
electrons  traveling  towards  the  electrode  are  reflected  by  the  moving  sheath  and  thus  gain 
energy  much  as  a  ping-pong  ball  would  after  being  struck  by  a  paddle.  Vendor  and 
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Boswell  provide  a  simple  calculation  based  on  their  numerical  data  and  a  hard  wall  model 
which  illustrates  the  degree  of  electron  energy  gain  in  electron-sheath  collisions.  In  this 
model  of  electron-sheath  collisions,  the  resulting  velocity  of  an  electron  incident  upon  the 
moving  sheath  is  given  by  (Vender  and  Boswell,  1990:728): 

^tlectron, final  “  ~^electron,milial'^'  ^  heath  (1-1) 

From  their  simulation  they  found  a  maximum  sheath  velocity  of  approximately  2  x  106 
m/s  and  the  velocity  of  a  thermal  electron  incident  upon  the  sheath  of  approximately  106 
m/s  (about  3  eV).  Using  equation  1.1  they  found  that  the  incident  electron  is  accelerated 
to  a  velocity  of  approximately  5  x  106  m/s.  This  corresponds  to  a  final  energy  of  roughly 
70  eV  (Vender  and  Boswell,  1990:728).  This  provides  a  simple  numerical  example  of 
the  role  of  the  sheath  and  electron  energy  gain  -  one  can  see  that  the  sheath  imparts  large 
amounts  of  energy  to  the  electrons.  This  energy  is  then  redistributed  within  the  discharge 
by  electron-neutral  collisions.  It  will  be  shown  below  that,  through  this  energy 
redistribution,  the  electrons  convert  "electric  power  into  chemical  power"  by  way  of 
ionization,  excitation,  and  dissociation  resulting  from  the  electron-neutral  collisions 
(Boeuf  and  Belenguer,  1990:  155).  This  energy  conversion  is  fundamental  to  many 
plasma  processing  applications,  but  before  moving  into  this  area  an  overview  of  the 
various  mechanisms  for  energy  deposition  within  the  discharge  is  warranted. 

The  electric  field  resulting  from  applied  RF  power  imparts  energy  to  charged  particles 
through  the  methods  outlined  above.  Since  the  electrons  are  lighter,  they  are  primarily 
responsible  for  the  redistribution  of  the  energy  within  the  discharge  (Gill  and  Donaldson, 
1931:723).  This  redistribution  takes  the  form  of  electron-neutral  collisions  which  may  be 
either  elastic  or  inelastic  (Flamm  and  Herb,  1989:23).  The  type  of  any  particular  collision 
depends  upon  the  gas  involved  and  the  electron  energy. 

In  an  elastic  collision,  the  sum  of  the  kinetic  energy  of  the  two  constituent  particles 
remains  a  constant  before  and  after  the  collision  (Goldstein,  1980: 1 17).  In  electron- 
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neutral  elastic  collisions  a  small  amount  of  energy  is  transferred  from  the  fast  electron  to 
the  slow  neutral  particle.  On  average,  this  energy  is  approximately  given  by  : 

(1-2) 

where  A  £  is  the  initial  energy  difference,  me  is  the  electron  mass,  and  Mn  is  the  neutral 
particle  mass  (McDaniel,  1964:24).  Since  me  «  Mn,  the  amount  of  energy  transferred  in 
any  individual  electron-neutral  elastic  collision  is  very  small.  Thus,  the  electron-neutral 
elastic  collision  is  an  inefficient  energy  transfer  mechanism  (Flamm  and  Herb,  1989:22). 
This  stated,  elastic  collisions  are  still  important  in  the  characteristics  of  the  discharge. 
First,  although  electron-neutral  elastic  collisions  are  inefficient,  the  energy  loss  suffered 
by  electrons  in  this  mechanism  may  be  significant,  since  they  can  occur  most  frequently 
in  a  discharge  (depending  upon  gas  type  and  operating  conditions).  Secondly,  elastic 
collisions  also  transfer  momentum.  This  tends  to  randomize  the  motion  of  the  electrons, 
further  increasing  their  importance  to  the  characteristics  of  the  discharge. 

In  an  inelastic  collision,  the  total  kinetic  energy  of  the  two  particles  is  not  conserved 
(Goldstein,  1980: 1 17).  In  electron-neutral  inelastic  collisions  electron  kinetic  energy 
goes  into  changing  the  internal  state  of  the  neutral  particle.  Since  this  energy  can  be  quite 
large,  this  is  often  an  efficient  mechanism  for  electron  energy  redistribution.  Electron- 
neutral  inelastic  collisions  describe  a  wide  class  of  collisions  which  include  ionization, 
excitation,  and  dissociation.  In  this  manner,  "electron  energy,  channeled  into  electron- 
neutral  inelastic  collisions,  maintains  the  supply  of  ions  and  radical  species"  in  the 
discharge  (Flamm  and  Herb,  1989:26).  The  ions  and  radical  species  are  the  chemical 
agents  used  in  many  plasma  processing  applications  (Flamm  and  Herb,  1989:26). 

In  electron-neutral  impact  ionization,  an  electron  provides  the  energy  required  to 
ionize  the  neutral  particle  by  removing  an  electron  from  the  neutral,  creating  an  electron 
and  a  positive  ion.  Any  electron  with  kinetic  energy  greater  than  the  ionization  potential 


5 


which  collides  with  a  neutral  gas  particle  has  a  chance  of  ionizing  that  neutral  particle. 
The  electron  loses  kinetic  energy  in  the  process.  The  electrons  that  have  gained  energy 
from  the  electric  field,  as  described  above,  provide  this  energy  to  ionization  in  order  to 
sustain  the  discharge,  creating  more  electrons  and  ions  which  are  continually  lost  to  the 
electrodes  and  walls  of  the  discharge.  As  it  will  be  shown  below,  the  ions  which  are  lost 
to  the  electrodes  are  important  in  many  plasma  processing  applications.  In  this  manner, 
the  electrons  provide  an  energy  conversion  and  transport  mechanism,  converting 
electrical  energy  from  the  field  into  the  energy  required  to  create  new  charged  particle 
species. 

New  particle  species  are  also  created  via  dissociation.  In  this  type  of  electron-neutral 
inelastic  collision  the  electron  provides  the  energy  required  to  produce  potentially 
chemically  reactive  species.  Again,  these  particles  are  important  in  many  plasma 
processing  applications,  providing  the  plasma  chemistry  to  perform  the  desired  process. 

In  an  electron-neutral  inelastic  collision  that  results  in  excitation,  the  neutral  particle 
is  excited  to  a  higher  energy  level.  The  electron  supplies  the  energy  required  to  excite  the 
neutral.  This  collision  process  is  important  in  reactors  where  photons  play  an  important 
role  in  the  either  plasma  processing  or  as  a  diagnostic  means  for  process  control,  since  the 
excited  state  neutrals  ultimately  decay  to  a  lower  energy  state,  thereby  emitting  a  photon 
in  the  process. 

Some  of  the  general  characteristics  of  the  RF  glow  discharge  have  been  reviewed. 
Next,  the  dependence  of  industry  on  plasma  processing  will  be  established,  as  well  as  the 
need  for  an  increase  in  the  scientific  understanding  of  these  devices  in  order  to  optimize 
their  use  and  aid  in  future  reactor  designs. 

Industrial  Plasma  Processing 

The  National  Research  Council  reports  that  the  following  industries  rely  heavily  on 
plasma  processing  of  one  sort  or  another  national  defense,  aerospace,  optics,  solar 
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energy,  telecommunications,  textile,  paper,  waste  management,  biomedicine,  automobile, 
and  microelectronics  industries  (National  Research  Council,  1991:9).  Shohet  emphasizes 
the  importance  of  plasma  processing  in  industry  by  estimating  the  dollar  amounts 
associated  with  the  potential  markets  for  some  of  the  most  important  industrial 
applications.  A  summary  of  his  estimates  are  shown  in  table  1.1  (Shohet,  1991:725). 


Application  Dollar  Value  ($  billion) 

Metal  Corrosion  Protection  50 

Plasma  Electronics  40 

Semiconductor  Processing  26 

High-Performance  Ceramics  5 

Tool  and  Die  Hardening  2 


Table  1.1.  Current  and  Potential  market  values  for  plasma  processing  (Shohet, 

1991:725). 

The  microelectronics  industry  is  perhaps  the  most  dependent  on  plasma  processing 
for  the  production  of  computer  chips.  One  specific  example  of  the  microelectronics 
industry's  dependence  on  plasma  assisted  processing  given  by  the  National  Research 
Council  is  the  construction  of  the  metallic  connections  on  an  integrated  circuit.  This 
process  is  comprised  of  the  following  steps  (all  of  which  depend  upon  plasma  processing 
to  some  degree):  (1)  deposition  of  a  layer  of  metal  on  a  substrate,  (2)  application  of  a 
pattern  on  the  metal  layer  (the  mask),  (3)  removal  of  the  metal  not  covered  by  the  mask 
leaving  the  desired  connections,  and  (4)  removal  of  the  mask.  Step  3,  the  etching 
process,  is  absolutely  dependent  on  plasma  processing  and  serves  well  to  illustrate  the 
need  for  an  increased  scientific  understanding  of  RF  glow  discharge  phenomena. 
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The  National  Research  Council  provides  a  contrast  between  chemical  etching  and 
plasma-assisted  etching.  Prior  to  the  advent  of  plasma-assisted  etching,  reactive 
chemicals  were  used  to  remove  the  material  not  protected  by  the  mask  (so  called  "wet" 
etching).  Unfortunately,  as  the  chemical  removes  the  unprotected  material,  the  material 
underneath  the  mask  is  exposed  to  the  chemical  and  is  also  etched  away  causing  the 
etched  pattern  to  take  on  the  appearance  in  figure  1.2.a.  This  feature  of  wet  etching  is 
highly  undesirable,  especially  as  the  size  of  integrated  circuits  decrease.  Only  plasma 
assisted  etching  (so  called  "dry"  etching)  is  capable  of  producing  the  high  fidelity  pattern 
transfer  required  for  today's  devices.  Figure  1.2.b  contrasts  the  capability  of  dry  etching 
over  wet  etching.  The  sharp  right  angle  etching  shown  in  figure  1.2.b  is  referred  to  as 
anisotropic  etching.  The  process  responsible  for  this  feature  will  be  explained  when  the 
etching  mechanisms  are  presented. 

Flamm  and  Herb  describe  an  early  parallel  plate  RF  glow  reactor  patented  by 
Reinberg  in  1975  (Flamm  and  Herb,  1989:65-66).  Figure  1.3  shows  a  schematic  of  the 
"Reinberg  style  reactor."  The  wafers  are  placed  on  the  lower  electrode,  which  may  be 
heated  or  cooled  to  enhance  the  chemistry  of  the  particular  plasma  processing  at  hand. 
This  system  was  originally  used  as  a  chemical  vapor  deposition  reactor,  and  is  now  one 
of  the  most  common  types  used  in  industry. 

In  his  doctoral  dissertation,  Barnes  describes  the  major  categories  of  plasma  assisted 
etching  processes  (Barnes,  1988: 14-16).  They  are  categorized  by  the  particles 
responsible  for  the  etching  process  and,  therefore,  provide  a  good  introduction  to  the 
basic  mechanisms  of  plasma  assisted  etching.  The  major  categories  are:  neutral  free 
radical  etching,  ion-assisted  etching,  and  sputter  etching  (Barnes,  1988:14-16). 

In  neutral  free  radical  etching,  free  radicals  formed  by  dissociation  in  the  discharge 
chemically  react  with  the  wafer.  This  chemical  reaction  forms  a  volatile  product  and 
results  in  an  etched  surface.  Note  that  this  is  similar  to  straight  chemical  etching  without 
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a  plasma  --  depending  upon  crystal  orientation,  there  may  not  be  a  preferential  direction 
of  the  etching  process  so  that  the  resulting  etch  pattern  resembles  those  found  in  figure 

1.2. a  (Barnes,  1988:14). 

Ion-assisted  etching  is  broken  into  two  subcategorizes:  ion-assisted  free  radical 
etching  and  chemically  enhanced  sputter  etching.  In  ion-assisted  free  radical  etching 
the  same  type  of  process  responsible  for  chemical  etching  is  at  work,  the  difference  being 
that  energetic  ions  supply  "activation  energy"  which  enhances  the  chemical  reaction 
(Sommerer,  1992:  2).  This  is  the  basis  for  the  anisotropic  etch  patterns  shown  in  figure 

1.2. b.  The  highly  anisotropic  ions  (with  velocities  primarily  directed  perpendicular  to  the 
electrode  surface  due  to  the  strong  axial  sheath  fields)  promote  etching  in  the  direction 
perpendicular  to  the  electrode  surface.  Anisotropic  etches  are  also  a  feature  of  chemically 
enhanced  sputter  etching ,  although  here  the  ions  perform  most  of  the  work  by  physically 
knocking  surface  atoms  off  the  wafer.  In  this  type  of  etching,  chemically  reactive  species 
tend  to  prepare  the  surface  (via  surface  heating  or  bond  breaking,  for  example)  so  as  to 
make  it  easier  for  the  ions  to  do  their  job.  This  reduces  the  requirement  for  high  energy 
ions  which  helps  reduce  the  potential  for  damage  to  the  wafer  (Barnes,  1988:15). 

Lastly,  straight  sputter  etching  employs  the  use  of  highly  energetic  (approximately 
500  eV)  ion  bombardment  to  dislodge  surface  atoms  from  the  wafer  (Barnes,  1988:15). 
This  type  of  etching  also  results  in  the  desirable  anisotropic  etching.  One  problem  with 
this  mechanism  is  that  the  high  energy  ions  can  result  in  surface  damage  to  the  wafer. 

Flamm  and  Herb  provide  a  diagram  (shown  in  figure  1.4)  which  summarizes  the 
plasma  assisted  etching  mechanisms  described  above.  Now  that  the  basic  mechanisms  of 
plasma  assisted  etching  have  been  outlined,  a  review  of  the  major  features  of  plasma 
etching  which  one  would  like  to  control  will  help  to  establish  the  need  for  an  increased 
scientific  understanding  of  the  discharge. 
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The  National  Research  Council  explains  that  in  any  plasma  etching  process,  one 
would  like  to  control  the  following  characteristics  (among  other  things):  anisotropy, 
selectivity ,  damage ,  and  uniformity  (National  Research  Council,  1991:21-23).  As 
mentioned  above,  anisotropy  refers  to  the  type  of  etching  shown  in  figure  1.2.b.  This 
type  of  etch  is  made  possible  by  the  presence  of  highly  directional  energetic  ions  which 
enhance  the  chemical  reaction  or  dislodge  surface  atoms.  Since  the  ions  are  accelerated 
through  the  high  sheath  field  toward  the  electrode,  controlling  the  structure  of  the  plasma 
sheath  is  a  prerequisite  for  controlling  ion  anisotropy.  Therefore,  scientists  would  like  to 
determine  the  sheath  structure  as  a  function  of  operating  parameters. 

Selectivity  refers  to  the  desire  to  remove  the  etch  layer  while  leaving  the  mask 
undamaged.  Scientists  recognize  the  fact  that  selectivity  depends  upon  both  ion  energy 
and  surface  chemistry  in  ion-assisted  etching.  However,  a  trade-off  exists  between 
selectivity,  anisotropy,  and  damage  —  if  the  same  ions  that  foster  anisotropy  are  too 
energetic  the  result  is  a  lack  of  selectivity,  not  to  mention  the  potential  for  damage  to  the 
wafer  due  to  the  high  energy  ions. 

Lastly,  uniformity  refers  to  the  ability  to  etch  over  a  wide  area  at  the  same  rate.  This 
ability  depends  upon  controlling  the  uniformity  of  ion  flux  (which,  it  has  been  pointed 
out,  depends  upon  the  sheath  structure).  This  may  be  further  complicated  by  localized 
effects.  For  instance,  the  specific  structure  of  an  individual  wafer  or  the  structure  of  a 
piece  of  a  single  wafer  may  change  the  etch  rate  from  wafer  to  wafer  or  across  a  single 
wafer. 

The  next  logical  question  is,  since  we  are  aware  of  some  of  the  characteristics  that 
influence  the  quality  of  plasma  etching  (and,  in  general,  any  plasma  assisted  process), 
how  can  these  characteristics  be  controlled?  Sommerer  describes  the  dilemma  plasma 
processing  practitioners  find  themselves  in  when  trying  to  influence  some  of  the  factors 
outlined  above: 
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The  reactor  conditions  may  be  optimized  for  the  desired  process  by  varying  the  fill 
gas  mixture,  flow,  pressure,  applied  voltage,  electrode  sizes,  and  substrate 
temperature.  Unfortunately,  typical  plasma  discharges  tend  to  respond  in  a  highly 
complicated  way  to  adjustments  of  any  one  of  these  parameters.  It  is  therefore 
difficult  to  adjust,  for  example,  the  ion  flux  to  the  wafer  or  substrate  without  also 
affecting  the  neutral  radical  flux.  This  complex  interaction  stymies  intuitive  attempts 
at  process  optimization,  and  motivates  much  of  the  effort  to  develop  more  rigorously- 
based  models  (Sommerer:  1992:2). 

Shohet  points  out  that  "the  key  to  further  advances  in  plasma-aided  manufacturing"  is 
the  understanding  of  the  complex  relationship  between  the  external  operating  parameters 
and  (in  the  case  of  plasma  assisted  etching)  the  characteristics  described  above  (Shohet, 
1991:  726).  To  date,  trial-and-error  has  been  the  method  used  when  attempting  to 
influence  many  plasma  processes  (Sommerer,  1992: 1).  In  light  of  the  important  role 
plasma  assisted  processing  plays  in  industry,  the  National  Research  Council  has 
recommended  a  plan  of  action  to  remedy  this  situation  (National  Research  Council, 

1991).  Paramount  in  this  plan  is  the  development  of  multi-dimensional  computer  models 
of  plasma  reactors  in  order  to  determine  the  relationships  between  external  operating 
parameters  and  plasma  processing  mechanisms  (National  Research  Council,  1991:51). 
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II.  RF  Glow  Discharge  Simulation 


General  Methods 

In  general,  there  are  two  primary  methods  for  simulating  the  RF  glow  discharge:  the 
fluid  method  and  the  particle  method.  This  section  briefly  discusses  the  major  features  of 
each.  Variations  of  these  basic  methods  exist,  as  well  as  methods  which  combine  features 
of  both.  These  techniques  are  not  discussed  here. 

Fluid  Method.  Boeuf  and  Belenguer  describe  the  fluid  method  which  they  have  used 
to  study  the  RF  glow  discharge.  In  the  fluid  method  the  modeler  begins  by  taking  a  finite 
number  of  successive  velocity  moments  of  the  Boltzmann  equation  in  order  to  find  the 
time  progression  of  the  discharge  particle  species.  Particle  species  movement  within  the 
discharge  is  accomplished  via  the  use  of  average  quantities  obtained  from  the  moment 
equations.  Assumptions  and  approximations  regarding  the  velocity  distribution  function 
of  each  species  must  be  made  in  order  to  close  the  system  of  equations  (Golant  and 
others,  1980:Ch  6). 

One  of  the  major  disadvantages  of  this  method  is  the  validity  of  making  apriori 
assumptions  regarding  the  particle  velocity  distribution  functions.  Another  short  coming 
of  this  approach  is  the  inability  to  predict  the  non-local  effects  important  in  RF  glow 
discharges.  Surendra  and  Graves  provide  an  example  (Surendra  and  Graves,  1991: 148). 

In  their  particle  simulation  of  a  one  dimensional  RF  discharge  in  helium  they  observed 
that  the  peak  electron  heating  rate  and  ionization  rate  did  not  occur  at  the  same  point  in 
the  discharge  due  to  the  long  mean  free  path  for  electron-neutral  collisions.  They  point 
out  that  fluid  models  are  hard  pressed  to  predict  this  type  of  non-local  behavior. 
Nevertheless,  fluid  model  results  similar  to  those  found  by  way  of  particle  based  methods 
have  been  obtained  in  much  shorter  computer  run-times  (Sommerer,  1992:8). 

Particle  Method.  In  the  particle  method  (sometimes  referred  to  as  the  particle-in-cell 
(PIC)  method  for  reasons  that  will  become  apparent)  the  trajectories  of  computer  particles 
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(called  super  particles)  representing  many  real  particles  are  simulated  (Sommerer, 

1992:4).  The  particle  method  makes  use  of  a  spatial  and  temporal  grid.  The  basic 
computational  cycle  is  shown  in  figure  2.1  (Birdsall,  1991:81).  In  this  method  Maxwell's 
equations  are  used  to  solve  for  the  fields  generated  by  the  charged  particles  and  the 
applied  voltage  on  the  electrodes.  Given  the  fields,  the  Newton-Lorentz  equation  of 
motion  is  then  used  to  update  the  trajectories  of  the  charged  particles  (Birdsall,  1991:65). 
As  noted  in  chapter  one,  charged  particle  collisions  are  important  in  RF  glow  discharges. 
Collisions  with  neutral  atoms  are  incorporated  via  the  Monte  Carlo  method.  The  fact  that 
the  particle  motion  and  calculation  of  the  fields  is  self-consistent  is  an  attractive  feature  of 
particle  simulation  —  it  "allows  the  model  to  'make  up  its  own  mind'  about  the  nature  of 
its  macroscopic  and  collective  behavior"  (Hockney,  1966: 1826).  Note  that  by  using 
Maxwell's  equations  and  a  spatial  grid,  the  requirement  of  calculating  the  force  on  each 
charged  particle  due  to  every  other  charged  particle  and  the  electrodes  is  negated  —  a 
computationally  impossible  task  for  even  a  moderate  number  of  particles  (Birdsall  and 
Langdon,  1991:9).  Furthermore,  in  simulating  the  discharge  we  are  interested  in  the 
collective  behavior  of  the  plasma  down  to  the  plasma's  characteristic  length,  the  Debye 
length.  By  using  the  spatial  grid  the  unwanted  detail  (not  to  mention  the  trouble 
associated  with  singularities  as  the  distance  between  particles  approaches  zero)  is  not 
calculated  (Birdsall  and  Langdon,  1991:9).  The  spatial  grid  is  used  to  calculate  the 
charge  and  current  densities  and  corresponding  fields  while  the  temporal  grid  is  used  to 
advance  the  simulation  from  time-step  to  time-step.  Both  grids  must  be  of  sufficient 
resolution  (small  enough  spacing  between  successive  grid  points)  to  maintain  numerical 
accuracy  in  the  simulation.  Furthermore,  the  spatial  grid  must  be  fine  enough  to  at  least 
resolve  a  Debye  length,  while  the  time-step  size  must  be  small  enough  to  resolve  the 
electron  plasma  frequency  (Hockney  and  Eastwood,  1988:32, 39, 333-334).  It  is 
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important  to  note  that  the  fields  are  "grid  quantities"  while  the  particle  positions  may  take 
on  continuous  values. 

Referring  to  figure  2.1,  the  particle  method  consists  of  five  major  calculations 
(Birdsall,  1991:81): 

(1)  Weight  the  particles  to  the  spatial  grid  to  form  the  densities; 

(2)  Calculate  the  fields  at  each  grid  point  given  the  densities  on  the  grid; 

(3)  Weight  the  fields  to  the  particle  locations  to  find  the  force  on  each  particle; 

(4)  Update  the  particle  positions  and  velocities  given  the  force  on  each  particle; 

(5)  Collide  charged  particles  with  neutrals. 

In  figure  2. 1  the  quantities  on  the  left-hand  side  of  the  arrow  are  transformed  by  the 
calculation  into  the  quantities  on  the  right-hand  side  of  the  arrow.  The  subscript  k  refers 
to  the  kth  particle,  while  the  subscript  j  refers  to  the  j*  grid  point.  Note  that  a  particle's 
position  is  not  updated  as  a  result  of  a  collision  --  only  its  velocity  is  changed.  The  series 
of  computations  is  carried  out  once  per  time-step. 

As  one  can  see  from  the  description  of  the  particle  method  its  major  advantage  is 
that  very  little  must  be  assumed  in  the  design  of  the  model  —  the  electric  fields,  particle 
production  and  loss,  and  particle  motion  are  calculated  using  first  principles.  As  alluded 
to  in  the  description  of  the  fluid  method,  particle  simulations  are  computationally 
intensive.  In  fact,  this  is  their  major  disadvantage.  This  disadvantage  is  becoming  less 
prohibitive  as  computer  speed  increases,  although  at  this  time  parameter  studies  of  two 
dimensional  particle  models  on  all  but  the  fastest  of  today's  computers  are  impractical. 
However,  this  does  not  preclude  the  need  for  two  dimensional  particle  models  of  RF 
discharges.  At  this  time,  a  two  dimensional  particle  model's  use  will  probably  be 
confined  to  a  limited  number  of  runs  in  order  to  "double  check"  fluid  model  results, 
validate  assumptions  and  approximations  made  in  the  fluid  models,  or  as  part  of  a  hybrid 
model  combining  features  of  both  methods.  Another  disadvantage  of  the  particle  method 
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is  numerical  heating.  In  order  to  reduce  the  computation  times  it  is  tempting  to  reduce 
the  number  of  simulation  particles  --  this  leads  to  artificial  heating  of  the  charged 
particles  due  to  noise  in  the  charge  density,  potential,  and  electric  field. 

Previous  Simulations 

To  date,  the  vast  majority  of  RF  glow  discharge  simulations  have  been  in  one  spatial 
dimension  (although  a  number  of  initial  papers  on  two-dimensional  models  were  included 
in  the  technical  program  of  the  45th  Annual  Gaseous  Electronics  Conference,  26  -  29 
October,  1992).  This  section  will  outline  a  few  of  them. 

Boeuf  and  Belenguer  simulated  one-dimensional  RF  glow  discharges  in  helium, 
nitrogen,  and  chlorine  using  the  fluid  method  (Boeuf  and  Belenguer,  1990:159).  They 
made  observations  regarding  the  electric  field,  charged  particle  densities,  current 
densities,  and  plasma  potential  at  various  points  in  an  RF  cycle  for  various  discharge 
conditions  and  geometries.  Special  emphasis  was  put  on  characterizing  the  mechanisms 
responsible  for  sustaining  the  discharge  (as  discussed  in  chapter  one)  for  various 
discharge  conditions.  Boeuf  and  Belenguer  used  a  "two  electron  group  fluid  model."  In 
this  method,  two  separate  groups  of  electrons  are  governed  by  their  own  set  of  (coupled) 
moment  equations.  One  group  represents  the  "beam  electrons"  which  gain  high  energy 
from  the  sheath  field  and  the  other  represents  the  slower  bulk  electrons  in  the  center  of 
the  discharge  (Boeuf  and  Belenguer,  1990: 158). 

Surendra  and  Graves  simulated  one-dimensional  parallel  plate  RF  glow  discharges  in 
helium  using  the  particle  method  with  Monte  Carlo  collisions  (Surendra  and  Graves, 
1991).  They  also  made  observations  regarding  the  discharge  sustaining  mechanisms 
(Surendra  and  Graves,  1991: 146-150).  A  result  of  their  simulation  runs  are  electron  and 
ion  energy  and  velocity  distribution  functions  for  various  discharge  conditions. 

Vender  and  Boswell  simulated  one-dimensional  parallel  plate  RF  glow  discharges  in 
an  atomic  hydrogen-like  gas  using  the  particle  method  with  Monte  Carlo  collisions 


15 


(Vender  and  Boswell,  1991).  They  found  that  the  scaling  of  the  sheath  length  and 
number  density  with  applied  voltage  was  consistent  with  the  Child-Langmuir  relation, 
namely 

3  1 

d,  «  V]  n}  (2.0) 

where  ds  is  the  sheath  length,  Vs  is  the  maximum  voltage  across  the  sheath,  and  nt  is  the 
electron  density.  In  their  simulation,  they  found  that  the  sheath  length  was  practically 

3 

independent  of  applied  RF  voltage,  so  that  ne  a  V2 .  They  also  obtained  electron  and  ion 
energy  distribution  functions  as  a  result  of  their  calculation.  This  work  was  based  on  a 
previous  model  by  Boswell  and  Morey  (Boswell  and  Morey,  1987). 

Unresolved  Problems 

There  are  many  unresolved  problems  concerning  RF  discharge  reactors  (in  terms  of 
their  application  to  plasma  processing)  which  have  not  yet  been  fully  addressed  by  either 
particle  or  fluid  models.  Some  of  these  problems  require  multidimensional  models  and 
are  thus  model  limited  since  multidimensional  RF  discharge  models  are  only  now 
beginning  to  be  developed.  Other  problems  are  data  limited  where  basic  cross-section 
and  gas  chemistry  data  is  required.  A  few  of  these  unresolved  problems  are  discussed 
below. 

Ion  Anisotropy  and  Uniformity.  The  characteristics  of  the  ion  flux  to  electrodes  is 
of  paramount  importance  since  they  play  a  pivotal  role  in  plasma  processing.  In  order  to 
examine  the  effects  of  bounding  surfaces  (container  walls)  on  the  uniformity  of  the  ion 
flux  across  an  electrode  surface  or  the  anisotropy  of  ions  bombarding  a  surface,  a  two 
dimensional  model  is  required.  Also  of  interest  is  the  ratio  of  radial  ion  flux  to  axial  ion 
flux.  This  requires  a  two-dimensional  model. 

Surface  Inhomogenieties.  Studying  the  effects  of  the  macroscopic  (wafer-electrode) 
and  microscopic  (within  an  individual  wafer)  surface  inhomogenieties  on  plasma 
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processing  characteristics  requires  at  least  a  two-dimensional  model.  Differences  in  the 
secondary  emission  coefficient  of  the  wafer  and  electrode  and  a  change  in  the  boundary 
conditions  on  the  potential  due  to  the  presence  of  the  wafer  have  not  yet  been 
investigated.  Likewise,  the  microstructure  of  the  wafer  effects  the  boundary  conditions 
on  the  potential,  thus  altering  the  fields  and  the  plasma  processing  characteristics  of  the 
discharge.  A  two  dimensional  model  is  required  to  study  these  effects. 

Gas  Chemistry  and  Surface-Gas  Chemistry.  In  general,  current  simulations  do  not 
model  the  neutral  gas  chemistry  (not  to  mention  the  more  complicated  gases  used  in  the 
plasma  processing)  nor  do  the  current  simulations  model  surface-gas  chemistry.  For 
instance,  much  is  unknown  about  the  role  of  excited  state  neutrals  in  the  discharge.  Some 
of  the  limitation  in  this  case  is  model  limited  and  requires  the  incorporation  of  neutral 
particle  diffusion  models  into  the  charged  particle  codes  and  some  may  be  data  limited 
for  some  gases  where  collision  cross-sections  may  be  unknown. 

Model  Development 

The  particle  method  with  Monte-Carlo  collisions  was  chosen  for  this  research. 
Although  the  inherently  long  computer  runs  are  unattractive,  a  particle  model  is  perhaps 
the  easiest  way  to  get  started  in  the  simulation  of  RF  discharges  due  to  the 
straightforward  modeling  techniques  and  the  lack  of  uncertainty  regarding  the  validity  of 
assumptions  and  approxinu*  Ions  required  in  the  fluid  method.  Ultimately,  a  suite  of 
models  consisting  of  particle  models  (single  and  multidimensional)  and  fluid  models 
(single  and  multidimensional)  is  an  ideal  goal.  Thus,  the  particle  models  can  be  used  to 
test  the  validity  of  the  assumptions  made  in  the  fluid  models.  On  the  other  hand,  the  fluid 
model  can  be  used  to  perform  parameter  studies  which  may  be  impossible  to  do  with  the 
particle  model.  Both  models  can  be  used  to  check  to  the  results  of  the  other.  Also, 
hybrid  particle-fluid  models  may  be  possible.  The  determination  of  how  to  partition  the 


model  can  be  based  on  either  particle  species  or  discharge  region  (sheath  vs.  bulk), 
depending  upon  which  is  more  suitable  to  either  method. 

Coordinate  System.  Given  the  decision  to  use  the  particle  method,  the  next  step  in 
designing  the  simulation  is  the  choice  of  coordinate  system.  In  order  to  address  topics 
such  as  those  outlined  in  chapter  one  (anisotropy,  uniformity,  etc.)  a  multidimensional 
model  is  required.  The  cylindrical  discharges  shown  in  figure  1.1  and  1.4  suggest  that 
cylindrical  coordinates  are  an  appropriate  choice.  Therefore,  two-dimensional  (2D) 
cylindrical  coordinates  (r-z,  no  theta  variations)  were  chosen.  This  coordinate  system  is 
attractive  in  two  respects.  First,  it  allows  one  to  study  the  radial  variations  across  the 
electrodes  (for  studying  anisotropy  and  uniformity)  while  also  being  able  to  compare  the 
axial  variations  to  existing  one-dimensional  models.  Secondly,  this  coordinate  system 
serves  as  a  good  starting  point  towards  a  full  three-dimensional  particle  simulation  of 
cylindrical  parallel  plate  discharges. 

The  particle  shape  in  a  particular  coordinate  system  is  conceptually  important  in  the 
development  of  the  model  equations.  As  an  example,  take  one-dimensional  cylindrical 
coordinates  in  r  only  (no  z  or  theta  variations).  The  particle  shapes  are  infinite  uniform 
cylindrical  shells,  centered  on  r  =  0  (Birdsall  and  Langdon,  1991:332).  Particle  motion 
corresponds  to  a  change  in  the  radius  of  this  infinite  cylinder  of  charge.  In  the  coordinate 
system  chosen  for  this  model,  the  particle  shapes  are  uniform  rings  of  charge  (Birdsall 
and  Langdon,  1991:333).  Particle  motion  consists  of  a  combination  of  axial  motion  and  a 
change  in  the  radius  of  the  ring  of  charge. 

Development  and  Discretization  of  the  Model  Equations  (Hockney  and  Eastwood, 
1988:26-32, 305-309).  The  next  step  in  the  development  of  the  model  is  the 
discretization  of  the  pertinent  equations.  Hockney  and  Eastwood  provide  a  method  for 
doing  so.  Once  the  applicable  set  of  equations  has  been  developed  the  particular 
algorithms  for  solving  the  equations  will  be  described.  As  stated  above.  Maxwell's 
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equations  are  used  to  solve  for  the  fields  and  the  Newton-Lorentz  equation  of  motion  is 
used  to  advance  the  velocity  and  position  of  the  particles  given  the  fields.  Maxwell's 
equations  are  solved  for  the  fields  given  the  charge  and  current  densities: 


div  B-0 

(2.1. a) 

,  -  -  1  dE 

c«r/fl-p0y+  , 

C  Ol 

(2.1b) 

div  E*^- 

(2.1.c) 

eo 

dB 

curl  E - 

(2.1.d) 

dt 

where  B  and  E  are  the  magnetic  and  electric  field  vectors,  respectively,  pis  the  charge 
density,  and  j  is  the  current  density.  The  force  on  a  particle  of  charge  q  and  mass  m 
under  the  influence  of  both  magnetic  and  electric  fields  is  given  by  the  Newton-Lorentz 
equation  of  motion: 

F-m--?(£+vxB)  (2.2) 

dt  v  ' 

where  v  is  the  particle's  velocity.  Equation  2.2  can  then  be  integrated  to  solve  for  the 
particle's  trajectory. 

In  the  electrostatic  plasma  approximation,  the  magnetic  field  due  to  the  moving 
charges  is  assumed  small  so  that  Maxwell's  equations  reduce  to  Poisson's  equation: 

VV-  (23) 

Co 

where  <j>  is  the  electrostatic  potential.  The  electric  field  is  then  given  by: 

E--grad<t>  (2.4) 

The  particle  trajectories  can  be  computed  by  integrating: 


19 


(2.5) 


dv  - 
m  —  -q  E 
dt  H 

and 


dx 

~dt 


-  v 


(2.6) 


Equations  2.3  through  2.6  are  the  model  equations  which  must  be  discretized  for 
solution  on  the  computer.  In  the  chosen  coordinate  system  Poisson's  equation  takes  on 
the  following  form: 


1  d  ( rd<t>\  _  ft?,  Z) 

r  dr\  dr)  dz 1  en 


(2.7) 


An  approximate  solution  to  equation  2.7  can  be  obtained  on  a  spatial  grid  with  uniform 
grid  spacing  Ar  and  Az  ( Ar  *  A z)  via  the  following  centered  finite  difference  equation 
(Swarztrauber  and  Sweet,  1975:45): 


+ 

(2.8) 

The  continuous  variables  of  equation  2.7  have  been  replaced  by  their  discrete 
counterparts.  The  potential  and  charge  density  are  now  grid  quantities  indexed  by  grid 
point.  The  index  /  refers  to  r  grid  points  and  the  index  j  refers  to  z  grid  points.  The 
grid  is  shown  in  figure  2.2.  Specific  software  for  solving  equation  2.8  is  discussed  below. 

The  right  hand  side  of  equation  2.8  is  formed  by  weighting  the  charge  from  the 
continuous  particle  positions  to  the  spatial  grid.  Next,  the  charge  weighted  to  the  gnd  is 
divided  by  the  appropriate  volume  element  to  form  the  charge  density  at  the  grid  points. 
The  notation  used  here  is  similar  to  that  used  by  Hockney  and  Eastwood  in  the 
development  of  a  one-dimensional  periodic  plasma  model  (Hockney  and  Eastwood, 
1988:3 1).  The  weighted  charge  at  each  grid  point  can  be  expressed  as: 
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(2.9) 


where  q  is  the  particle's  charge,  Ns  is  the  number  of  "real"  charged  particles  per 
computer  (or  super)  particle,  Np  is  the  total  number  of  particles,  and  w[rk  tzk,rt,Zj )  is 

some  function  which  weights  the  charge  at  the  particle  positions  ( rk,zk )  to  the  grid 
positions  (r,,zy).  In  practice,  the  sum  over  all  particles  in  equation  2.9  is  generally  done 
only  over  those  particles  in  the  grid  cells  centered  on  the  grid  point  (i,  j).  Specific 

methods  for  doing  the  charge  weighting  are  discussed  below.  The  divisor  in  equation  2.9 
is  the  volume  of  a  grid  cell  centered  on  grid  point  (4  j)  (Birdsall  and  Langdon,  1990:333- 

334).  In  2D  r-z  cylindrical  coordinates  the  volumes  are  rings  Az  by  A r .  Referring  to 
figure  2.3,  for  i  *  0  the  volume  elements  are  given  by: 


and  for  i  -  0 , 


V,, -2*/-(Ar)(Az) 


V0J-^(Ar)J(Az) 


(2.10) 


(2.11) 


Once  the  potential  is  found  on  the  grid,  the  electric  field  components  can  be  found  by 
finite  differencing  equation  2.4  which,  in  our  coordinate  system,  takes  on  the  form: 


F  * 

E-  T 

Et  - 

dz 


(2.12) 


Equations  2.12  have  the  following  centered  finite  difference  approximations  (Hockney 
and  Eastwood,  1988:  30): 


F  m  6-1  ■/  ~ 

%  2(Ar) 

F  m  6.7-1  ~ 

^  m  2  (A z) 


(2.13) 
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Now  that  the  field  equations  have  been  put  in  their  discrete  forms,  a  method  for 
integrating  the  equations  of  motion  must  be  decided  upon  in  order  to  discretize  the 
equations  of  motion  (equations  2.5  and  2.6).  The  leapfrog  finite  difference 
approximation  offers  an  efficient  way  to  integrate  equations  2.5  and  2.6  (Hockney  and 
Eastwood,  1988:28).  It  is  called  the  leapfrog  method  because  the  position  and  velocity 
are  known  at  alternating  half  time-steps  as  shown  in  figure  2.4  (Birdsall  and  Langdon, 
1991: 13).  The  positions  and  forces  on  the  particles  are  known  at  integer  time-steps  while 
the  particle  velocities  are  know  at  half  time-steps.  The  particle  positions  and  velocities 
are  updated  via  (Hockney  and  Eastwood,  1988:32): 


F'J*) 

-  -  vr  +  — 12 

*  Ns  rr^ 


-  v. 


z*  '*/2  N,  mk 


v,+vk(a'> 


(2.14) 


where  k  is  the  particle  index,  A /  is  the  size  of  the  time-step,  t  is  the  time-step  index,  and 
F  is  the  force  on  the  particle.  Note  that  since  the  force  is  known  on  the  grid,  it  must  be 
weighted  back  to  the  particle  (with  the  same  weighting  used  for  forming  the  charge 
density). 

Normalization  of  the  Model  Equations  (Hockney  and  Eastwood,  1988:33  -  34). 
Now  that  the  model  equations  have  been  discretized,  some  simplifications  can  be  made 
by  the  introduction  of  dimensionless  and  normalized  variables.  Starting  with  the  leapfrog 
equations  of  motion,  if  one  introduces  the  following  dimensionless  variables  (tilded 
quantities  are  either  dimensionless  or  normalized): 
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*  Fr  (AO"  (Z15) 

a'm  N,m(Ar) 

a--FM 
1  Nsm(Az) 

where  ar  and  a:  are  dimensionless  acceleration  components,  the  equations  of  motions 
simplify  to  (Hockney  and  Eastwood,  1988:33): 


v,  -  v,  +  aT 

\,+y2  \,-y2 

v,  -  v  +  a, 

lk.,*y2  lk.,-y2 


rk.,+ l-h'+Vr 


k.t+y 

k,i*y2 


(2.16) 


Note  that  integrating  the  equations  of  motion  has  been  reduced  to  four  additions  per  time- 
step  (aside  from  the  particle  weighting). 

In  order  to  complete  the  process,  the  electric  field  components  will  be  normalized. 
First,  substituting  the  expression  for  weighting  the  electric  field  components  to  the 
particles  (the  force  on  the  particle)  into  the  normalized  acceleration  components  gives  the 
following: 
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mk  (A r) 

q‘‘\'2w{r"z*'r"zj)Ezl,,t  •(A/): 


rk,t 


a,  - 

ZM 


(2.17) 


mk  (Az) 

Equations  2.13  and  2.17  suggest  the  following  normalization  for  the  electric  field 
components: 

\2 


p  (aQ  /  _  \ 

"2m,  (Ar)2 

f  g* ML fi  \ 

“  2  m,  (Ar)2  *'J+U  > 


(2.18) 


2  mk  (Az) 

Now  that  the  model  equations  have  been  simplified  it  is  useful  to  summarize  the  final 
equations  and  their  use  before  moving  on.  First,  equation  2.9  through  2.1 1  are  used  to 
form  the  charge  density  on  the  grid  (using  specific  weighting  schemes  discussed  below). 
Second,  Poisson's  equation  (equation  2.8)  is  solved  to  find  the  potential  on  the  grid. 
Equations  2. 18  are  then  used  to  find  the  normalized  electric  field  components  on  the  grid. 
The  leapfrog  integration  method  is  then  used  to  update  each  particle's  position  and 
velocity  components.  Rewriting  equations  2. 16  in  terms  of  the  normalized  electric  field 
components  gives: 


(2.19) 


rk.,.X-rk,+Vr 


I-Z*.,+Vz 


In  practice,  the  sum  over  all  grid  points  in  equation  2. 19  is  done  over  only  the  nearest  one 
or  four  grid  points,  depending  upon  the  weighting  method  used.  This  sequence  is 
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accomplished  once  per  time-step.  The  above  equations  do  not  completely  describe  the 
model.  The  solution  of  Poisson's  equation,  boundary  conditions,  particle  and  field 
weighting,  and  charged-neutral  particle  collisions  must  be  discussed. 

Solution  of  Poisson's  Equation  on  the  Grid.  Efficient  FORTRAN  subprograms  for 
solving  equation  2.8  (among  other  separable  elliptic  partial  differential  equations)  have 
been  developed  by  Swarztrauber,  Sweet,  and  Adams  and  made  available  by  the  National 
Center  for  Atmospheric  Research  via  the  internetwork  in  the  FTSHPAK  software  package 
(Adams  and  others,  1988).  The  F1SHPAK  subprograms  are  used  to  solve  equation  2.8  on 
the  grid.  FISHPAK  uses  a  cyclic  reduction  algorithm  to  solve  the  block  tridiagonal 
system  of  equations  that  result  from  equation  2.8.  The  general  method  is  described  in 
(Sweet,  1977)  while  the  actual  computer  subprograms  are  documented  in  (Swarztrauber 
and  Sweet,  1975). 

The  FISHPAK  subprogram  HWSCYL  was  used  to  solve  equation  2.8.  (Swarztrauber 
and  Sweet,  1975:37  -  59).  This  subprogram  requires  that  the  right-hand  side  of  equation 
2.8  reside  on  the  grid  prior  to  calling  the  subprogram.  When  the  subprogram  completes 
the  grid  contains  the  solution  to  equation  2.8  in  terms  of  the  potential.  In  addition,  the 
subprogram  HWSCYL  requires  that  boundary  conditions  be  specified.  The  boundary 
conditions  can  be  of  the  following  form:  solution  specified,  normal  derivative  of  the 
solution  specified,  or  unspecified  at  r  =  0  (due  to  symmetry).  In  this  case  three  boundary 
conditions  (at  r  =  R,  z  =  0,  and  z  =  L)  are  required,  where  R  is  the  radius  of  the  discharge 
tube  and  L  is  the  distance  between  the  electrodes.  The  boundary  at  r  =  0  is  allowed  to 
remain  unspecified  since,  by  symmetry,  the  normal  component  of  the  derivative  of  the 
solution  vanishes. 

Boundary  Conditions .  The  grid  used  in  the  calculation  and  the  boundaries  upon 
which  conditions  must  be  specified  are  shown  in  figure  2.5.  Since  the  voltage  on  each 
electrode  is  known  at  all  times,  the  solution  can  be  specified  for  the  boundary  conditions 
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at  z  =  0  and  z  =  L.  Either  or  both  electrodes  may  be  held  at  a  constant  voltage  or  driven 
at  some  frequency.  The  voltages  on  the  electrodes  are  specified  by: 


l'«  -  *  K,0  Slrl[-  .1  f„a  1  (A»)] 

Vl  -  Kl  *  V„L  si"[2  *  f„L  t  (A/)] 


(2.20) 


where  V0  and  VL  are  the  voltages  on  the  electrodes  at  z  =  0  and  L  respectively,  and 
VdCL  are  dc  bias  voltages,  Vr/o  and  V rf^  are  the  peak  driving  voltages,  and  fT/^  and  fr/^ 

are  the  driving  frequencies  in  Hz. 

The  boundary  condition  at  the  glass  wall  is  more  complicated.  If  we  treat  the  glass 

wall  as  "an  infinite  sticky  insulator"  after  Hockney  we  can  specify  the  normal  derivative 

of  the  solution  at  the  remaining  boundary  (Hockney,  1966: 1829).  Any  charged  particle 

which  strikes  the  glass  wall  sticks  to  the  wall  --  its  charge  is  weighted  to  the  grid  point(s) 

adjacent  to  the  cell  and  deleted  from  the  simulation,  leaving  behind  a  surface  charge.  It 

can  be  shown  (see  appendix  A)  that  the  normal  component  of  the  electric  field  within  a 

semi-infinite  dielectric  with  infinite  dielectric  constant  vanishes.  This  will  serve  as  the 

boundary  condition  imposed  behind  the  surface  charge  which  forms  on  the  glass  wall. 

dd> 

Hence,  the  normal  derivative  of  the  solution  will  be  specified:  Er  - - -  0.  Figure  2.6 

dr 


depicts  the  grid,  glass  wall,  surface  charge,  and  boundary  condition. 

In  addition  to  specifying  the  conditions  on  the  potential  at  the  boundaries,  the  particle 
behavior  at  the  boundaries  must  also  be  specified.  The  particle  behavior  at  r  =  R  has 
already  been  specified  above.  If  a  particle  crosses  the  grid  line  at  r  =  0,  by  symmetry,  the 
r  components  of  the  particle's  position  and  velocity  are  negated.  The  z  components  are 
left  unchanged.  At  the  electrodes  the  particle  behavior  is  specified  by  reflection  and 
emission  coefficients  after  Surendra  and  Graves  (Surendra  and  Graves,  1991: 145).  A 
charged  particle  which  comes  into  contact  with  one  of  the  electrodes  may  be  absorbed  or 
reflected  by  the  wall  depending  upon  the  comparison  of  a  uniform  random  number  to  the 


specified  reflection  coefficient.  Particles  that  are  not  absorbed  by  the  wall  undergo 
specular  reflection  back  into  the  discharge.  Furthermore,  a  charged  particle  which  strikes 
the  wall  may  cause  the  emission  of  an  electron  from  the  wall.  Again,  a  secondary 
electron  is  created  depending  upon  the  comparison  of  a  uniform  random  number  to  the 
specified  emission  coefficient.  Emitted  electrons  are  given  an  initial  velocity  directed 
perpendicular  to  the  wall  with  some  user  specified  energy.  In  general,  ions  are  given  a 
zero  reflection  coefficient  while  electrons  are  given  a  zero  emission  coefficient. 

One  last  word  concerning  the  boundaries.  In  the  non-boundary  regions  of  the  grid  the 
electric  field  components  are  found  by  simply  applying  equation  2. 18  to  all  interior  grid 
points.  On  the  radial  boundaries  the  electric  field  is  specified  by  the  boundary  conditions 
themselves  (as  described  above).  On  the  electrode  boundaries  applying  equation  2.18 
presents  a  problem  in  that  the  grid  "runs  out"  and  equation  2.18  can  not  be  used.  A 
simple  solution  used  by  Vender  and  Boswell  is  to  estimate  the  normal  component  of  the 
electric  field  at  the  boundary  by  "linear  extrapolation  from  the  two  values  nearest  the  grid 
point"  (Vender  and  Boswell,  1990:726).  This  method  is  adopted  for  use  here.  Of  course, 
the  component  of  the  electric  field  parallel  to  the  electrode  at  the  axial  boundaries  is  zero, 
since  the  electrode  is  an  equipotential  surface. 

Particle  and  Field  Weighting.  Birdsall  and  Langdon  describe  several  methods  for 
charge  and  field  weighting  (Birdsall  and  Langdon,  1991: 19  -  20, 308  -  309, 336  -  337). 
The  charge  weighting  schemes  will  be  described  here  with  the  understanding  that  the 
same  method  is  used  for  weighting  the  fields  back  to  the  particles.  The  goal  is  to  weight 
the  charge  from  the  continuous  particle  locations  to  the  discrete  grid  points  to  form  the 
right  hand  side  of  Poisson's  equation.  Any  of  three  weighting  schemes  may  be  used: 
nearest- grid-point  (NGP),  bilinear,  or  volume  weighting  (bilinear  in  r2).  Birdsall  and 
Langdon  refer  to  the  bilinear  in  r2  method  as  area  weighting.  Since  the  weighting  is 
actually  based  on  a  cylindrical  volume,  the  name  volume  weighting  seems  more 
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appropriate  and  will  be  used  here,  although  the  algorithm  is  similar  to  that  described  in 
Birdsall  and  Langdon  (Birdsall  and  Langdon,  1991:336  -  337). 

NGP  weighting  assigns  all  of  a  particle's  charge  to  the  nearest  grid  point.  This 
method  is  efficient,  since  it  only  requires  the  particle's  rounded  position  components. 
However,  the  use  of  NGP  results  in  a  noisy  plasma  —  as  a  particle  moves  across  a  cell  the 
charge  assigned  to  a  single  mesh  point  suddenly  jumps  from  zero  to  the  particle's  charge 
(Birdsall  and  Langdon,  1991:20). 

Figure  2.7  will  help  to  explain  both  the  bilinear  and  volume  weighting  schemes 
(Birdsall  and  Langdon,  1991:309).  In  bilinear  and  volume  weighting  the  particle's 
charge  is  assigned  to  the  four  nearest  grid  points  (Birdsall  and  Langdon,  199i:336). 
Referring  to  figure  2.7,  the  proportion  of  the  charge  assigned  to  grid  point  A  is  (Birdsall 
and  Langdon,  1991:336): 

_ Area  or  Volume  of  region  a _  (2  21) 

Area  or  Volume  of  regions  (a  +  b  +  c  +  d) 

The  remaining  three  grid  points  are  handled  in  a  similar  manner. 

Monte  Carlo  Collisions.  Charged  particle-neutral  collisions  are  accomplished  via 
Monte  Carlo  methods  with  the  aid  of  the  null  collision  technique.  Birdsall  provides  a 
general  overview  of  incorporating  Monte  Carlo  collisions  into  particle  simulations 
(Birdsall, 1991:79-82).  Birdsall  points  out  that  (unlike  the  development  so  far)  the 
computer  particles  are  treated  as  individual  particles  in  the  collision  process.  So  far,  the 
computer  particles  have  been  treated  as  many  physical  particles  (Birdsall,  1991:80). 

In  the  null  collision  technique,  a  "fake"  total  collision  cross-section  is  constructed 
which  must  conform  to  the  following  rules:  1)  it  must  be  greater  than  the  sum  total  of  all 
other  collision  cross-sections  for  all  energies  and  2)  it  must  be  proportional  to  the  inverse 
of  the  cl.arged  particle's  speed.  In  this  way,  the  "fake"  total  collision  frequency  is 
independent  of  particle  energy  and  constant  for  the  entire  simulation  so  that  the 
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calculation  of  a  cross-section  for  each  particle  at  each  time-step  is  precluded  (Birdsall, 
1991:82).  Now  that  the  fake  total  collision  cross-section  has  been  constructed,  the  null 
collision  technique  can  be  carried  out. 

The  first  step  is  to  calculate  the  percentage  of  particles  to  collide  each  time-step: 

p-  1  -  exp[-«<ai  V  ( A/)]  (2.22) 

where  ngas  is  the  neutral  gas  number  density,  is  the  fake  total  collision  cross-section, 
and  v  is  the  particle  speed  (Birdsall,  1991:82).  Since  the  fake  total  cross-section  is  given 
by: 


where  k  is  a  proportionality  constant,  equation  2.22  can  be  calculated  independent  of 
particle  speed  and  at  the  start  of  the  simulation.  Given  the  percentage  of  the  total  number 
of  particles  to  collide  with  neutrals,  a  random  subset  is  constructed  from  the  total  set  of 
charged  particles.  The  percentage  is  calculated  at  the  start  of  the  simulation,  while  the 
random  subset  is  determined  each  time-step.  The  subset  of  particles  to  collide  is  selected 
as  follows.  A  "skipping  constant"  is  calculated  by  taking  the  inverse  of  equation  2.22 
(this  is  calculated  at  the  start  of  the  simulation)  (Smith,  1991).  A  uniform  random 
number  Ri  (0  <,  Rj  £  1)  is  drawn  and  multiplied  by  the  skipping  constant  at  the  start  of 
each  time-step.  This  gives  a  new  index  of  the  first  particle  in  the  list  of  charged  particles 
to  collide.  Subsequent  particles  are  chosen  by  looping  through  the  list,  incrementing  the 
loop  index  by  the  skipping  constant  each  time  through  the  loop.  This  procedure  produces 
a  random  subset  of  the  total  set  of  charged  particles. 

Now  that  the  subset  of  particles  to  collide  is  known,  the  type  of  collision  must  be 
determined.  This  is  done  by  comparing  the  collision  cross-sections  relative  to  the  fake 
total  cross  section  to  a  uniform  random  number,  R2  (0  £  R2  £  1)  as  in  Hockney  and 
Eastwood  (Hockney  and  Eastwood,  1988:379).  Wherever  R2  falls  in  the  sequence  of 
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relative  collision  cross-sections  is  the  collision  process  that  is  carried  out.  An  example 
will  help  to  illustrate  the  method.  Say,  for  instance,  one  would  like  to  model  two  real 
collisions.  Including  the  null  collisions  this  makes  a  total  of  three  possible  collision 
processes.  First,  calculate  the  following  quantities: 


J  jaU 


m. 


ai  +  cr2 


* fake 


(2.24) 


m3  -  l-(m,  +  wj 

Note  that  the  mi  sum  to  one,  with  Q<mi<m1<mi<\.  Next,  draw  the  uniform  random 


number  R2.  Where  R2  falls  in  the  sequence  of  quantities  calculated  via  equations  2.24 
determines  the  collision  which  occurs.  For  OsR,  <  m, ,  real  scattering  process  one 
occurs.  For  m1  s  R,  <  m^,  real  scattering  process  two  occurs.  If  m2sR,  si,  then  a  null 
collision  occurs  -  nothing  is  done  (Hockney  and  Eastwood,  1988:379). 

Once  we  have  calculated  which  particles  undergo  which  collision  processes  during 
each  time-step,  the  only  detail  in  the  collision  process  left  to  discuss  is  the  collision 
kinetics.  Before  describing  the  collision  kinetics  for  the  collision  processes  included  in 
the  model,  the  number  of  velocity  components  carried  in  the  simulation  must  be 
mentioned.  As  described  above,  the  model  equations  of  motion  require  only  two  velocity 
components  --  one  for  each  dimension  simulated.  However,  in  order  to  properly  account 
for  the  particle's  kinetic  energy,  the  third  velocity  component,  ve ,  is  also  carried.  The 
quantity  vg  is  only  altered  by  collision  processes.  This  is  similar  to  the  technique  used  by 
Burger  in  his  simulation  of  a  one-dimensional  plasma  diode  (Burger,  1967:660-661).  As 
mentioned  above,  the  particle  positions  are  not  affected  by  a  collision.  Only  the  particle 
velocities  are  perturbed. 

The  present  simulation  includes  electron-neutral  elastic,  electron-neutral  ionization 
and  ion-neutral  charge  exchange  collisions.  The  incorporation  of  the  remaining  collisions 
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is  straight  forward,  given  the  basic  framework  constructed  for  the  existing  collisions.  For 
elastic  collisions,  the  scattered  electrons  energy  is  calculated  via  (McDaniel,  1964:24): 

E  ■  f  * 

electron,  scattered  electron  incident 

In  the  case  of  ionization,  the  ionization  potential  for  the  particular  gas  being  simulated  is 
subtracted  from  the  total  energy  of  the  incident  electron.  This  energy  is  then  randomly 
distributed  between  the  scattered  electron  and  the  electron  created  by  ionization.  The 
created  ion  is  given  the  energy  of  a  thermal  ion.  For  an  ion-neutral  charge  exchange 
collision,  the  ions  energy  is  reduced  to  the  user  defined  thermal  ion  energy.  Once  the 
energy  of  the  constituent  particles  is  calculated,  it  is  distributed  to  the  particle's  three 
velocity  components  isotropically.  That  is,  all  directions  of  motion  are  made  equally 
probable.  Although  the  collision  kinetics  have  been  greatly  simplified,  they  provide  a 
start  towards  a  more  realistic  representation.  Specific  suggestions  for  improvement  will 
be  made  in  chapter  four. 

A  discussion  of  the  method  for  specifying  the  collision  cross-sections  completes  the 

section  on  particle  collisions.  Birdsall  describes  a  method  for  specifying  electron-neutral 
collision  cross-sections  (Birdsall,  1991:81).  The  quantities  ,  EL 

ionization  ionization 

E  ,  and  E,  in  figure  2.8  are  used  to  describe  an  approximation  to  the 
electron-neutral  ionization  cross-section.  Similarly,  the  quantities  o ^  ^  and  E^  ^ 

in  figure  2.9  are  used  to  describe  an  approximation  to  the  electron-neutral  elastic  collision 
cross-section.  In  addition,  the  fake  total  collision  cross-section  is  described  by  the 
proportionality  constant  in  equation  2.23,  chosen  to  meet  the  criteria  listed  above  for  the 
fake  total  collision  cross-section.  The  proportionality  constant  should  be  chosen  to 
minimize  the  number  of  null  collisions  by  making  the  fake  total  collision  cross-section  as 
close  as  possible  to  (but  always  greater  than)  the  real  total  collision  cross-section. 


(2.25) 

I 
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Although  not  used  here,  another  method  used  by  Ramos  for  increasing  the  efficiency 
of  the  null  collision  technique  is  the  use  of  multiple  fake  collision  cross  sections  in  order 
to  minimize  the  number  of  null  collisions  over  different  energy  regions  (Ramos,  1990:55- 
58).  As  more  collision  processes  are  added  to  the  simulation  it  may  be  worth  the  added 
time  savings  to  incorporate  this  method. 

Initial  Particle  Loading.  The  are  many  choices  for  the  method  of  loading  the  initial 
positions  and  velocities  of  the  simulated  electrons  and  ions.  One  method  suggested  by 
Birdsall  is  to  load  the  particles  uniformly  in  position  and  use  a  normalized  thermal 
velocity  distribution  (Birdsall,  1991:390-391).  This  method  was  attempted  with 
unfavorable  results.  At  first  blush  it  seems  that  almost  any  reasonable  (devoid  of  large 
oscillations)  initial  phase  space  portrait  would  suffice  (Vender  and  Boswell,  1990:726). 
Unfortunately,  due  to  current  memory  and  simulation  time  constraints,  the  simulation 
must  be  run  with  a  limited  number  of  particles  (<  1.6  million  per  species  on  the  Silicon 
Graphics  IRIS  workstation  used  for  this  research).  This  limitation  effects  the  allowable 
initial  positions.  If  the  particles  are  placed  within  the  simulation  uniformly  in  position,  a 
large  number  of  the  electrons  (which  are  much  faster  than  the  relatively  immobile  ions) 
initially  flee  to  the  bounding  surfaces.  Given  enough  electrons  and  simulation  time,  this 
would  not  present  a  problem.  Regrettably,  due  to  the  small  number  of  particles,  a  large 
percentage  of  the  total  number  of  electrons  are  lost  and  it  is  unclear  whether  the 
simulation  could  ever  "straighten  itself  out"  In  order  to  offset  this  problem  the  particles 
are  loaded  so  that  fewer  particles  are  put  in  the  sheath  regions,  adjacent  to  the  electrodes 
and  glass  wall.  Furthermore,  a  number  of  electrons  may  be  initially  "stuck"  to  the  glass 
wall  (typically  a  few  tens  of  electrons  per  grid  point)  to  impede  initial  electron  loss  to  the 
glass  wall. 
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The  initial  particle  velocities  are  loaded  using  a  normalized  thermal  velocity 
distribution  as  in  Birdsall  (Birdsall,  1991:390-391).  The  initial  speed  for  each  particle  is 
given  by: 

V  “  V*ern»l  V~2  l°8  (^)  (2.26) 

where  vtkermal  is  a  user  defined  thermal  speed  and  ^  is  a  uniform  random  number. 

Special  care  must  be  taken  to  exclude  any  -  0 .  The  speed  calculated  via  equation  2.26 
is  then  divided  isotropically  between  the  velocity  components  (Birdsall,  1991:390-391). 

Computer  Code.  The  computer  program  which  implements  the  model  described 
above  was  written  in  the  C  programming  language.  The  program  runs  on  SUN 
SparcStations  and  Silicon  Graphics  workstations.  The  software  does  not  take  advantage 
of  any  system  specific  features  so  that  porting  the  program  to  other  computer  models 
should  be  straightforward. 

The  source  code  consists  of  five  C  source  files  each  providing  a  major  piece  of 
functionality.  In  addition,  each  source  file  has  its  own  header  file  (excluding  rf.c)  that 
contains  an  interface  to  the  functions  provided  by  that  source  file  (as  well  as  any  required 
constants,  special  data  types,  etc.).  Each  of  the  functions  contained  in  the  source  files  has 
a  short  description  of  its  purpose,  input  parameters  and  output  parameters.  The  five 
packages  are:  particles. c,  grid.c,  input _data.c,  collisions. c,  and  rf.c.  The  source  code  file 
rf.c  is  the  main  program,  while  the  four  other  packages  provide  the  functionality  their 
names  suggest  In  addition,  a  header  file  named  pic.h  provides  general  constants  used  by 
all  of  the  packages. 

partlcles.c.  This  package  defines  data  types,  allocates  storage  for,  and  provides 
functions  for  manipulating  the  particle  lists,  loading  the  initial  particle  positions  and 
velocities,  and  moving  the  particles  each  time-step. 

input_data.c.  This  package  provides  a  function  for  reading  the  input  data  as  well  as 
global  variable  definitions  for  the  input  data.  A  function  for  calculating  some  constants 
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based  on  the  input  data  is  also  provided.  All  user  input  is  done  via  file.  The  format  of  the 
input  file  is  described  below. 

gricLc.  This  package  defines  data  types,  allocates  storage  for,  and  provides  functions 
for  operating  on  the  simulation  grid.  Functions  for  charge  weighting,  solving  Poisson's 
equation  (an  interface  to  the  FISHPAK  routine),  and  calculating  the  electric  field  is 
provided.  Since  FISHPAK  was  written  in  FORTRAN  and  the  simulation  is  written  in  C, 
special  care  was  taken  in  writing  the  function  that  interfaces  to  FISHPAK.  Since 
FORTRAN  stores  two-dimensional  arrays  in  column  major  format  and  C  stores  them  in 
row-major  format,  the  grid  must  be  converted  to  and  from  column  major  format  each 
time-step. 

collisions.*:.  This  package  provides  the  functionality  to  perform  Monte  Carlo 
collisions.  Functions  are  provided  to  calculate  the  collision  cross  sections  as  a  function  of 
particle  energy.  A  function  is  provided  for  determining  collision  types  and  collision 
kinetics  based  on  the  cross  section  functions  and  particle  energy. 

Input  File  Format.  A  sample  input  file  is  shown  in  appendix  C.  Table  2.1  indicates 
the  meaning  of  each  variable  name  as  well  as  the  units  associated  with  these  quantities. 


Input  Variable 

Explanation  (Units) 

dz 

z  grid  cell  size  (m). 

Nz 

Number  of  grid  cells  in  z.  Nz  *  dz  gives  distances  between 
electrodes  (none). 

dr 

r  grid  cell  size  (m). 

Nr 

Number  of  grid  cell  in  r.  Nr  *  dr  gives  radius  of  discharge 
tube  and  electrodes. 

dt 

Time-step  size  (seconds). 

Nt 

Number  of  time  steps  to  simulate  (none). 

npr 

npz 

The  quantities  npr  and  npz  define  the  maximum  number  of 
particles  per  cell  =  npr  *  npz. 
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MiMe 

Ve 

Vi 

P 

Ns 

sigma_total 


sigma_max_ion 

eOion 

elion 

e2ion 

sigma_max_elas 

eOelas 

sigma_ch_ex 


SecE 


Seel 


VSec 

RefE 

Refl 

VrfO 

frfO 

VdeO 

VrfL 

frfL 

VdcL 

weighting 


Nout 

Pout 

Gout 


The  ratio  of  ion  mass  to  electron  mass  (none). 

Speed  of  a  thermal  electron  (m/s). 

Speed  of  a  thermal  electron  (m/s). 

Neutral  gas  pressure  (Torr). 

Number  of  real  particles  per  super  particle  (none). 

The  proportionality  constant  in  equation  2.23  used  to 
describe  the  fake  total  collision  cross  section. 

The  next  five  parameters  describe  the  electron-neutral 
ionization  cross  section  as  in  figure  2.8. 


The  next  two  parameters  describe  the  electron-neutral 
elastic  cross  section  as  in  figure  2.9. 

Constant  ion-neutral  charge  exchange  collision  cross- 
section  (mA2). 

Coefficient  to  describe  secondary  electron  emission  due 
electron  bombardment.  Always  set  to  0  (none). 

Coefficient  to  describe  secondary  electron  emission  due 
ion  bombardment  (none). 

Velocity  of  a  secondary  electron  (m). 

Electron  electrode  reflection  coefficient  (none). 

Ion  electrode  reflection  coefficient  (none). 

The  next  six  parameters  are  the  constants  in  equation  2.20 
used  to  describe  the  voltages  on  the  electrodes 
(voltages  in  volts  and  frequencies  in  Hz). 


Weighting  scheme  to  use.  1  =  NGP.  2  =  Bilinear. 

3  =  Volume.  Other  =  Error,  (none). 

Diagnostic  output  line  every  Nout  time-steps  (none). 

Write  potential  on  grid  to  file  every  Pout  time-steps  (none). 

Write  electron  phase  space  file  every  Bout  time-steps 
(none). 
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lout  Write  ion  phase  space  file  every  Eout  time-steps  (none). 

Fout  Write  electric  field  components  on  grid  to  file  every  Fout 

time-steps  (none). 

Sout  Write  charge  on  glass  wall  to  file  every  Sout  time-steps 

(none). 

Cout  Write  charge  on  grid  to  file  every  Cout  time-steps  (none). 

Rout  Write  the  right-hand-side  of  Poisson's  equation  to  file  every 

Rout  time-steps  (none). 

Tout  Update  ionization  diagnostic  file  every  Tout  time-steps 

(none). 

seed  Random  number  generator  seed.  Change  seed  will  start  a 

new  pseudo-random  number  sequence  (none). 

Table  2.1.  Input  File  Description. 

Program  Compilation,  Linking  and  Execution.  The  program  is  compiled  and 
linked  (through  a  command  found  in  the  source  code  directory)  by  typing  build  at  the 
operating  system  command  line.  This  produces  an  executable  file  called  rf.  In  order  to 
execute  the  program  the  following  command  line  is  executed: 

rf  input_file  [electron_phasc_space_flle]  [ion_phase_space_file]  [wall_charge_file]  <return> 
where  input_file  is  of  the  format  specified  above.  The  command  line  arguments  in  square 
brackets  are  optional.  They  specify  phase  space  files  and  wall  charge  from  previous  runs. 
This  feature  may  be  used  to  "warm  start"  the  simulation  from  where  a  previous  run 
stopped. 

Output  Files.  The  program  produces  a  number  of  output  files.  Each  output  file 
name  has  the  following  format:  input_file.time_step.output_fi!e.  The  input_file 
corresponds  to  the  input  file  name  specified  on  the  command  line.  The  time-step 
corresponds  to  the  time-step  at  which  the  output  file  was  created.  The  output_file 
corresponds  to  which  output  file  this  is.  The  following  quantities  may  be  written  to  file 
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during  a  simulation:  (1)  potential  at  each  grid  point,  (2)  electron  phase  space,  (3)  ion 
phase  space,  (4)  electric  field  components  at  each  grid  point,  (5)  surface  charge  on  the 
glass  wall,  (6)  weighted  charge  at  each  grid  point,  and  (7)  the  right  hand  side  of  Poisson's 
equation  at  each  grid  point.  In  addition  a  file  is  created  and  updated  which  stores  the 
location  of  ionizing  collisions  that  have  occurred  during  the  simulation.  Small  programs 
were  written  to  analyze  the  data  contained  in  the  phase  space  file  to  obtain  the  average 
particle  velocity  components,  average  particle  energy,  and  average  number  of  particles 
per  cell.  These  programs  will  not  be  described. 
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III.  Model  Verification,  Experimental  Results  and  Discussion 


Single  Particle  Test 

Two  simple  cases  involving  the  judicious  placement  of  a  single  particle  within  the 
simulated  discharge  tube  were  run  to  test  the  simulation  (Smith,  1991).  In  case  1,  a  single 
electron  was  placed  on  the  radial  axis  of  the  discharge  tube  and  exactly  midway  between 
the  two  electrodes.  In  the  second  case,  a  single  electron  was  again  placed  on  the  radial 
axis  of  the  tube,  but  displaced  to  one  side  of  the  midway  point  between  the  electrodes.  In 
both  cases  the  particle  was  given  zero  initial  velocity  and  both  electrodes  were  held  at 
zero  potential.  Both  cases  have  an  analytic  solution  for  the  motion  of  the  electron  and  are 
thus  good  tests. 

Case  1.  A  single  electron  with  zero  initial  velocity  placed  on  the  radial  axis  of  the 
discharge  tube  and  midway  between  two  grounded  electrodes  will  remain  at  rest.  This 
can  be  shown  via  the  method  of  images  (Feynman  and  others,  1964:  6-9  to  6-10).  The 
electron  has  two  image  charges  +e°  and  +e°2  located  as  shown  in  figure  3. 1.  These 
image  charges  also  have  images.  The  image  charge  +e°  has  an  image  -e\.  Similarly, 
the  image  charge  as  an  image  -e\ .  These  images  of  image  charges  go  on  forever. 
The  motion  of  the  electron  follows  from  the  force  on  the  electron  due  to  the  infinite 
number  of  image  charges. 

The  force  exerted  on  a  charged  particle  by  a  second  charged  particle  is  given  by 
Coulomb's  law: 


p  .  - 

rn  •» 


4  n  s0r 


1  r 


(3.1) 


where  Fl2  is  the  force  on  particle  1  due  to  particle  2,  g,  is  the  charge  on  particle  1,  Q2  is 
the  charge  on  particle  2,  e0  is  the  permittivity  of  free  space,  r  is  the  distance  between  the 
charges,  and  er  is  a  unit  vector  pointing  from  Q2  to  Q  •  In  the  case  at  hand,  the  force 
exerted  on  the  electron  due  to  the  infinite  number  of  image  charges  is  zero,  since  the 


38 


force  due  to  each  pair  of  image  charges  vanishes  due  to  the  symmetrical  location  of  the 
images  about  the  electron:  The  force  on  the  electron  due  to  q[  exactly  cancels  the  force 
on  the  electron  due  to  q'2,  where  q  takes  on  values  of  +e  or  -e.  Therefore,  it  has  been 
shown  that  a  single  electron  with  zero  initial  velocity  placed  on  the  radial  axis  of  the 
discharge  tube  and  midway  between  two  grounded  electrodes  will  remain  at  rest. 


The  simulation  was  run  with  the  conditions  described  above  for  50  time  steps.  As 
expected,  the  particle  did  not  move  --  it  remained  on  axis  and  midway  between  the  two 
electrodes  for  the  duration  of  the  simulation.  This  is  a  good  qualitative  test  of  the  field 
solver  and  particle  mover.  Case  2  provides  a  more  quantitative  test  and  insight  into  some 
of  the  model's  limitations. 

Case  2.  The  simulation  was  run  with  the  same  conditions  as  case  1  except  the 
electron  was  moved  slightly  off  the  midway  point  between  the  two  electrodes.  The 
acceleration  of  the  electron  as  a  function  of  position  may  be  calculated  analytically  (again 
using  the  method  of  images)  and  compared  to  the  approximate  value  calculated  and  used 
by  the  simulation. 

Figure  3.2  shows  the  position  of  the  electron  along  with  six  of  the  infinite  number  of 
image  charges.  The  electron  is  a  distance  z  from  the  left  electrode  and  consequently  a 
distance  L  -  z  from  the  right  electrode.  We  would  like  to  find  an  analytic  expression  for 
the  acceleration  of  the  electron  as  a  function  of  the  coordinate  z.  Let  the  image  charges 
q'„  be  labeled  in  the  same  fashion  as  case  1.  Let  r\  be  the  distance  from  image  charge  q\ 
to  the  electron  located  at  z.  The  force  on  the  electron  due  to  image  charge  q'„  is  given  by: 
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The  total  force  on  the  electron  due  to  the  infinite  number  of  images  is  given  by: 
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The  distance  between  each  image  charge  and  the  electron,  r'n ,  can  be  expressed  in  terms 

of  the  location  of  the  electron,  z,  and  the  distance  between  the  electrodes,  L,  and  is  shown 
beneath  each  image  charge  in  figure  3.2.  Upon  substitution  of  the  first  twelve  r'n  into 

equation  3.3  and  3.4  and  simplification: 


-  e2  r  1  1  1  1 

"l6we0[  z2+(L-zY~(L  +  z)2+(2L-z) 


(3.4) 


Since  er  lies  along  the  z-axis  for  every  image  charge,  the  acceleration  of  the  electron  due 
to  all  the  image  charges  is  entirely  along  the  z-axis.  Upon  simplification  of  equation  3.4 
and  the  use  of  Newton's  second  law,  the  z-component  of  the  acceleration  can  be  written 
as: 


Equation  3.5  provides  an  exact  analytic  solution  for  the  acceleration  of  the  electron  which 
may  be  used  to  test  the  accuracy  of  the  approximate  acceleration  calculated  by  the 
simulation. 


The  acceleration  used  by  the  simulation  to  advance  the  particle  from  time-step  to 
time-step  may  be  found  from  the  model  output  and  equations  of  motion.  The  z- 
component  of  the  electron's  velocity  is  updated  each  time  step  via: 


where  v^l+V2  is  the  z-component  of  the  particle's  velocity  at  time-step  t  +  1/2,  is 
the  z-component  of  the  particle's  velocity  at  time-step  t  -  \l2,az  l  is  the  z-component  of 


the  particle's  acceleration  at  time-step  t,  and  A t  is  the  size  of  the  time-step.  Equation  3.6 
can  be  rearranged  to  find  the  z-component  of  the  particle's  acceleration  at  time-step  t: 
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(3.7) 


"  A/ 

To  compare  the  analytic  solution  with  that  of  the  model,  the  simulation  was  run  for  35 
time-steps.  The  particle's  position  and  velocity  components  were  appended  to  an  output 
file  after  each  time-step.  The  analytic  solution  was  calculated  at  each  time-step  using  the 
model  calculated  position  of  the  electron  and  thirty  terms  in  the  infinite  summation  of 
equation  3.5.  Note  that  the  analytic  solution  to  the  acceleration  is  not  used  to  separately 
update  the  particle's  position.  In  this  way,  the  only  errors  considered  are  single  time-step 
errors  in  the  acceleration  --  cumulative  effects  in  the  trajectory  of  the  electron  due  to 
errors  in  the  calculation  of  the  acceleration  at  each  time-step  are  not  calculated.  Also  note 
that  each  "electron"  in  the  simulation  is  actually  modeled  as  10^  real  electrons.  That  is, 
each  model  electron  is  scaled  to  have  a  charge  of  10^  e  and  a  mass  of  10^  me  (where  e  is 
the  charge  of  an  electron  and  me  is  the  mass  of  an  electron).  The  acceleration  used  by  the 
model  was  calculated  at  each  time-step  using  the  model  calculated  velocity  and  equation 
3.7. 

Figure  3.3  shows  the  percentage  difference  between  the  analytically  calculated 
acceleration  and  the  model  calculated  acceleration.  It  is  interesting  to  note  that  the 
percent  difference  reaches  local  maxima  and  minima  as  the  electron  travels  towards  the 
electrode  (in  this  case  the  initial  displacement  was  towards  the  left  electrode  so  that  the 
electron  has  an  acceleration  in  the  negative  z  direction).  Figure  3.4,  a  plot  of  the 
electron's  normalized  z-component  of  position  versus  time-step,  provides  a  clue  not  only 
to  the  overall  source  of  error,  but  also  to  the  source  of  the  maxima  and  minima  in  the 
percent  difference  between  the  two  solutions.  Note  that  spatial  grid  points  exist  at  integer 
values  of  normalized  position. 

From  figure  3.3,  the  first  minima  in  the  percent  difference  occurs  at  time-step  ten. 

One  can  see  from  figure  3.4  that  this  corresponds  to  the  electron  being  closest  to  a  grid 
point.  The  same  is  true  of  time-step  twenty,  and  so  on.  As  the  electron  approaches  a  grid 


41 


point  the  charge  that  is  weighted  from  the  electron  to  the  grid  and  the  acceleration  that  is 
weighted  from  the  grid  to  the  electron  is  most  accurate.  Although  the  percent  difference 
improves  near  a  grid  point,  why  is  the  error  still  quite  large?  The  answer  is  in  the  bilinear 
weighting  scheme  used  in  this  case.  Although  bilinear  weighting  (as  discussed  in  chapter 
two)  reduces  the  noise  associated  with  particle  and  force  weighting,  it  is  ill-suited  for  the 
problem  at  hand.  When  the  charge  due  to  the  electron  is  weighted  to  the  two  nearest  grid 
points  the  problem  is  changed  from  a  single  particle  to  two  particles  (whose  individual 
mass  and  charge  sum  to  that  of  the  original  electron)  located  one  grid  space  apart. 
Although  this  "smearing"  effect  is  permissible  (and  desirable)  when  dealing  with  the 
collective  behavior  of  a  plasma,  it  introduces  errors  (as  seen  in  figure  3.3)  into  the 
calculation  when  dealing  with  a  single  particle. 

One  way  to  get  around  this  two-for-one  particle  problem  is  by  using  the  NGP 
weighting  scheme.  Since  the  particle  is  maintained  as  a  single  particle  whose  charge  is 
weighted  wholly  to  the  nearest  grid  point,  one  would  expect  the  difference  between  the 
analytic  solution  and  the  model  solution  to  decrease  as  the  electron  approaches  a  grid 
point.  The  results  for  NGP  weighting  with  otherwise  the  same  conditions  as  above  are 
shown  in  figures  3.5  and  3.6.  Note  that  when  the  particle  is  close  to  a  grid  point  the 
percent  difference  between  the  analytic  and  model  solution  improves  over  bilinear 
weighting.  However ,  when  the  particle  gets  far  from  a  grid  point  the  error  grows  to 
greater  than  thirty  percent  —  much  worse  than  bilinear.  This  is  the  result  of  the  whole  of 
a  particle's  charge  being  weighted  to  a  grid  cell  which  can  be  up  to  half  a  grid  cell  from 
the  actual  position  of  the  particle. 

The  results  of  both  single  particle  tests  give  confidence  in  the  proper  functioning  of 
the  model.  Furthermore,  case  2  provides  some  insight  into  a  limitation  of  the  model.  In 
particular,  the  simulation's  strength  is  not  in  modeling  single  particle  effects. 
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Experimental  Results  and  Discussion 

Noise  Problem.  A  problem  with  the  two-dimensional  model  prevents  running  the 
model  to  convergence  and  obtaining  quantitatively  meaningful  data.  This  problem  is 
associated  with  noise  in  the  potential  in  the  grid  cells  on  radial  axis.  This  section  explains 
the  cause  of  this  noise.  Also,  a  method  is  developed  which  allows  quantitatively 
meaningful  analysis  of  some  of  the  model  results. 

Before  attributing  the  cause  of  the  noise,  a  description  of  the  problem  is  presented. 
The  circumstances  surrounding  the  noise  on  the  grid  are  as  follows.  During  the  first  ten 
RF  periods  some  noise  in  the  potential  is  seen  on  the  radial  axis  of  the  grid.  This  noise 
does  not  extend  very  far  towards  the  radial  wall  (perhaps  a  few  grid  cells).  This  initially 
small  amount  of  noise  in  the  potential  corresponds  to  a  much  larger  than  normal  electric 
field  in  the  bulk.  This  gives  rise  to  an  anomalous  heating  of  the  electrons  in  the  bulk. 
During  the  course  of  the  simulation,  this  heating  artificially  increases  the  ionization  rate 
in  the  bulk  on  radial  axis.  This  rise  in  plasma  density  tends  to  increase  the  noise,  which 
in  turn  further  increases  the  plasma  density.  By  the  beginning  of  the  eleventh  RF  cycle 
the  noise  in  the  potential  is  quite  large  and  the  plasma  density  on  axis  is  growing  by  leaps 
and  bounds.  Figure  3.7.a  shows  a  plot  of  the  potential  as  a  function  of  position  within  the 
discharge  tube  at  the  beginning  of  the  eleventh  RF  cycle.  Figure  3.7.b  shows  a  one¬ 
dimensional  slice  of  the  same  potential  on  the  radial  axis.  By  this  time  the  on-axis 
potential  is  quite  noisy. 

As  explained  above,  the  model  uses  a  uniform  radial  grid  to  compute  the  field 
quantities.  This  is  believed  to  be  the  cause  of  the  noise  and  can  be  explained  as  follows. 
The  volume  elements  used  to  calculate  the  charge  density  at  grid  points  on  the  radial  axis 
are  eight  times  smaller  than  the  adjacent  volume  elements  just  one  radial  grid  cell  away 
(see  equation  2. 10  and  2. 1 1).  Since  the  number  of  simulation  particles  in  the  two  grid 
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cells  are  approximately  equal,  this  disparity  in  the  size  of  the  volume  elements  leads  to  a 
disparity  in  the  charge  densities  associated  with  the  two  grid  points  --  the  charge  density 
of  an  on-axis  grid  point  is  eight  times  larger  than  that  found  one  radial  grid  cell  away. 
This  large  increase  in  charge  density  is  believed  to  give  rise  to  the  noise  in  the  potential. 
This  claim  is  supported  by  figure  3.8  which  shows  four  plots  of  potential  as  a  function  of 
position  which  were  generated  by  loading  sixteen  cells  with  equal  numbers  of  electrons 
and  ions.  The  particles  were  loaded  randomly  in  position  within  each  grid  cell 
(therefore,  each  grid  cell  is  overall  neutral).  The  electrodes  are  held  at  zero  potential.  In 
each  plot  the  same  number  of  real  particles  in  a  grid  cell  is  maintained  by  varying  the 
number  of  simulation  particles  and  the  number  of  real  particles  per  simulation  particle. 
For  instance,  if  the  number  of  simulation  particles  is  increased  by  a  factor  of  ten  the 
number  of  real  particles  per  simulation  particle  is  decreased  by  a  factor  of  ten.  Figures 
3.8.a  -  3.8.c  shows  the  effect  of  increasing  the  number  of  computer  particles  on  the  noise 
on  the  grid  for  sixteen  on-axis  grid  cells.  As  expected,  as  the  number  of  computer 
particles  is  increased  the  noise  in  the  potential  is  reduced.  Figure  3.8.d  should  now  be 
compared  to  figure  3.8.c.  For  equal  numbers  of  computer  particles  the  on-axis  noise  in 
the  potential  is  a  factor  of  ten  greater  than  the  off-axis  noise.  This  data  supports  the 
above  explanation  for  the  source  of  the  noise. 

Given  the  above,  running  the  present  two-dimensional  model  to  convergence  is 
impossible.  Furthermore,  the  results  for  times  prior  to  the  large  increase  in  number 
density  and  potential  (in  addition  to  not  being  converged)  are  quantitatively  unusable  due 
to  the  anomalous  heating.  Nevertheless,  after  the  simulation  has  run  for  ten  RF  cycles 
many  of  the  qualitative  features  of  the  discharge  are  evident  and  can  be  commented  on. 

In  order  to  do  some  quantitative  analysis  (short  of  incorporating  a  non-uniform  radial  grid 
or  going  to  a  cartesian  geometry)  the  model  was  converted  to  a  one-dimensionaJ  axial 
model.  This  was  accomplished  in  short  order  by  using  only  one  component  of  position 
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(axial  distance  between  the  two  electrodes)  and  two  velocity  components  (an  axial 
component  and  a  component  perpendicular  to  the  axial  comoonent).  The  solution  of 
Poisson's  equation  is  simplified  in  one  dimension.  Most  other  components  of  the 
simulation  were  used  without  (or  with  only  minor)  modification.  In  this  respect, 
dropping  back  to  one-dimension  for  quantitative  analysis  provides  a  good  method  for 
validating  most  of  the  code  given  the  present  problems  of  running  in  two  dimensions. 
Also,  the  one-dimensional  model  represents  a  discharge  were  the  radial  boundary  is 
assumed  to  be  located  at  r  =  »,  thus  providing  a  limiting  solution  for  comparison  to  the 
two-dimensional  model  with  a  finite  radial  boundary.  After  all,  one  of  the  primary 
motivations  in  moving  to  a  two-dimensional  model  is  to  determine  the  effect  of  the  radial 
boundary  on  a  discharge's  plasma  processing  characteristics  (via  ion  motion  in  the 
sheaths).  As  the  discharge  tube  radius  is  reduced,  one  would  expect  its  effect  on  the 
sheath  ions  to  increase.  Once  the  noise  in  the  two-dimensional  model  is  corrected,  the 
effects  of  varying  the  radius  of  the  discharge  can  be  determined. 

The  remainder  of  this  chapter  is  organized  as  follows.  First,  the  experimental 
conditions  for  the  simulation  runs  are  described.  Second,  one-dimensional  results  are 
presented  with  comparison  to  previous  models  and  theory.  Third,  two-dimensional 
results  will  be  shown  caveated  by  the  above  discussion.  Quantitative  analysis  is 
performed  on  the  one-dimensional  data  since  it  does  not  suffer  from  the  noise  problem. 
The  qualitative  features  of  the  two-dimensional  data  are  highlighted  with  comparisons  to 
the  one-dimensional  runs.  Furthermore,  comment  is  made  on  ion  motion  in  the  sheaths 
as  it  pertains  to  plasma  processing  applications  (ion  anisotropy  and  uniformity). 

Although  no  conclusions  can  be  drawn  from  this  specific  two-dimensional  data,  it 
illustrates  the  relevant  issues  and  motivates  further  development  of  the  model. 

Experimental  Conditions.  The  discharge  discussed  in  this  section  (unless  otherwise 
stated)  consists  of  a  helium-like  gas  within  a  discharge  tube  with  4  cm  electrode 
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separation.  For  two-dimensional  runs  the  discharge  tube  radius  is  2  cm.  The  ion  nr  ass  is 
7344  me  and  the  collision  cross-sections  closely  approximate  those  of  helium.  The 
electron-neutral  collision  cross-sections  are  based  on  experimental  data  compiled  by 
Kieffer  (Kieffer,  1973: 1-7).  A  constant  (independent  of  energy)  ion-neutral  charge 
exchange  collision  cross-section  of  1.5  x  10-19  m2  is  used  (Surendra  and  Graves, 

1991: 145).  The  left  electrode  (z  =  0)  is  driven  sinusoidally  at  lOMhz  between  +/-500  V 
with  no  DC  bias  (see  equation  2.20).  The  right  electrode  (z  =  4  cm)  is  held  at  zero 
potential.  The  neutral  gas  pressure  is  250  mtorr.  The  electrons  and  ions  are  loaded  as 
discussed  in  chapter  two  (ions  with  a  thermal  speed  of  1.094  x  103  m/sec,  electrons  with  a 
thermal  speed  of  1.186  x  106  m/sec).  For  ions  the  electrodes  are  perfectly  absorbing  — 
when  an  ion  strikes  an  electrode  it  is  always  absorbed.  On  the  other  hand,  electrons  are 
reflected  in  25%  of  their  collisions  with  the  electrodes.  Charged  particle  bombardment  of 
the  electrodes  causes  no  secondary  particle  emission. 

The  non-physical  parameters  (grid  and  time-step  sizes)  were  chosen  to  satisfy  the 
following  relations  (Hockney  and  Eastwood,  1988:32): 


wp  (A/)<<  2 
A  S  An 


(3.8) 


where  wp  is  the  electron  plasma  frequency,  A t  is  the  time-step  size,  A  is  the  grid 
spacing,  and  AD  is  the  Debye  length.  The  discharges  considered  here  have  a  plasma 
density  of  a  few  10 15  nr3  and  bulk  electron  temperatures  of  a  few  eV.  This  gives  a 
Debye  length  of  a  few  10-4  m  and  an  electron  plasma  frequency  of  a  few  109  rad/sec . 
For  the  one-dimensional  runs  a  grid  spacing  of  104  m  and  time-step  0.1  nsec  were 
chosen  to  satisfy  relations  3.8  (due  to  memory  and  run-time  constraints  a  grid  spacing 
four  times  larger  than  this  was  used  in  the  two-dimensional  runs).  At  ten  MHz  driving 
frequency  the  choice  of  time-step  size  conveniently  gives  1000  time-steps  per  RF  cycle. 
In  order  to  maintain  the  discharge  geometry  400  grid  cells  are  used.  In  the  one- 
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dimensional  runs  anywhere  from  10,000  -  40,000  simulation  particles  per  species  were 
used.  Although  the  data  was  less  noisy  when  more  simulation  particles  were  used  the 
physical  quantities  (potential,  plasma  density,  etc.)  were  comparable  (for  the  same 
conditions)  over  all  ranges  of  number  of  simulation  particles  used. 

The  one-dimensional  simulations  were  run  for  250  -  500  RF  cycles.  Tracking  the 
total  number  of  simulation  particles  over  time  gave  an  initial  indication  of  convergence  — 
when  the  total  number  of  particles  remained  constant  over  a  couple  of  hundred  RF  cycles 
(within  one  or  two  percent  of  the  total  number  of  particles)  the  simulation  was  assumed 
converged.  In  the  case  of  the  one-dimensional  model  the  limiting  factor  in  attaining 
convergence  is  the  time  it  takes  for  the  ions  in  the  center  of  the  discharge  to  reach  the 
electrodes  (Vender  and  Boswell,  1990:725).  Simulations  started  with  a  plasma  density 
which  is  too  low  or  too  high  demand  longer  run  times  dictated  by  the  rate  of  loss  or 
production.  One  can  use  the  time  it  takes  for  an  ion  in  the  center  of  the  discharge  to  reach 
an  electrode  as  an  estimate  of  a  time  characteristic  of  convergence.  For  the  conditions 
stated  above  this  corresponds  to  approximately  5.0  psec  (50  RF  cycles).  Thus,  a 
characteristic  time  for  reaching  convergence  is  some  number  of  ion  transit  times  (-*50  RF 
cycles)  dictated  primarily  by  how  close  the  initial  bulk  plasma  density  is  to  the  converged 
solution. 

One-Dimensional  Results.  Two  one-dimensional  simulation  cases  are  presented. 
First,  the  results  of  a  one-dimensional  run  with  the  parameters  described  above  are 
discussed.  Second,  the  results  are  presented  from  a  one-dimensional  run  with  the  same 
parameters  as  the  first,  except  that  no  ion-neutral  collisions  are  included.  A  comparison 
of  the  two  runs  provides  a  good  example  of  the  importance  of  ion  collisions  on  the  sheath 
properties  of  the  RF  glow  discharge. 

Figures  3.9  and  3. 10  show  the  potential  and  electric  field  (respectively)  as  a  function 
of  position  within  the  discharge  at  various  times  in  the  RF  cycle.  Note  the  action  of  the 
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sheaths  in  response  to  the  electrode's  changing  potential.  The  sheath  at  the  left  electrode 
vanishes  at  tot  /  2n  =  0.25  while  the  sheath  at  the  right  electrode  vanishes  at  cut  /  2x  = 

0.75,  where  to  is  the  RF  frequency  in  rad/sec,  and  t  is  time  in  sec. 

Figure  3.11  shows  the  charged  particle  densities  as  a  function  of  position  within  the 
discharge  at  three  times  in  the  RF  cycle.  As  previously  discussed,  the  electrons  are 
excluded  from  the  region  adjacent  to  an  electrode  at  those  times  in  the  cycle  when  a 
sheath  is  formed  between  the  electrode  and  the  plasma.  As  the  sheath  contracts  the 
electrons  enter  the  previously  forbidden  region.  When  the  sheath  has  completely 
vanished  some  electrons  reach  the  electrode  as  indicated  in  the  non-zero  electron  number 
density  at  the  left  electrode  in  figure  3. 1  l.c.  The  average  electron  energy  as  a  function  of 
position  within  the  discharge  for  various  times  in  the  RF  cycle  is  shown  in  figures  3.12. 
Much  can  be  determined  from  figure  3.12  about  the  response  of  the  electrons  to  the  time- 
varying  field  in  the  discharge.  Figure  3. 12.a  shows  the  average  electron  energy  at  the 
beginning  of  the  RF  cycle.  Note  that  the  energy  is  not  symmetric  -  the  average  electron 
energy  near  the  right  plasma-sheath  boundary  is  greater  than  the  energy  near  the  left 
plasma-sheath  boundary.  This  can  be  explained  by  following  the  electron  energy  through 
a  half-cycle.  Figures  3. 12.b  through  3. 12.e  show  the  average  electron  energy  at  cot  /  2  ji  = 
0. 125, 0.25, 0.375,  and  0.5.  It  can  be  seen  that  at  cot  /  2n  =  0.5  the  same  lack  of  symmetry 
is  apparent,  except  that  the  higher  energy  electrons  are  near  the  left  electrode.  Let's 
follow  the  electrons  into  in  the  sheath,  starting  with  figure  3. 12.a.  At  this  time  the 
electrons  which  have  made  it  into  the  region  near  the  left  electrode  have  (on  average) 
higher  energies  than  the  bulk  electrons.  At  first  thought  this  seems  unreasonable.  After 
all,  (referring  to  figure  3.10.a)  there  is  no  field  in  the  sheath  region  at  this  time  that  will 
support  accelerating  the  electrons  to  such  a  high  energy.  It  turns  out  that  only  those  high 
energy  elections  (previously  in  the  bulk)  make  it  into  the  left  most  portion  of  the 
discharge  at  this  time  in  the  cycle.  At  cut  /  2ji  =  0.25  the  slower  electrons  have  had  a 


chance  to  enter  the  region  and  the  average  energy  decreases  slightly.  At  this  point  the 
sheath  is  about  to  reform.  As  the  sheath  advances  (figures  3.12.d  and  e)  the  electrons  in 
the  sheath  are  swept  out  of  the  sheath  region,  gaining  energy  from  the  sheath.  This 
explains  the  lack  of  symmetry  at  the  beginning  and  middle  of  the  RF  cycles.  Between  the 
middle  and  the  beginning  of  the  cycle  the  entire  process  is  repeated,  except  at  the  right 
electrode.  Consequently,  those  electrons  at  the  right  plasma-sheath  boundary  are  left  with 
higher  average  energies  (figure  3.12.a).  This  finding  (and  explanation)  as  well  as  an 
observed  average  bulk  electron  energy  of  approximately  5  eV  and  plasma  density  in  the 
center  of  the  discharge  of  3.4  x  10 15  m  3  is  consistent  with  the  work  of  others  under 
similar  conditions  (Surendra  and  Graves,  1991:146,  150). 

Figure  3. 1 1  also  shows  that  the  ion  density  is  not  modulated  by  instantaneous  changes 
in  the  electric  field.  Even  as  the  sheath  vanishes  and  reforms  the  ion  density  remains 
fairly  constant.  This  is  due  to  the  large  mass  of  the  ions  as  compared  with  that  of  the 
electrons  -  at  this  high  driving  frequency  the  ion  density  responds  to  a  period  average 
field  which  continually  draws  ions  to  the  electrodes. 

Figure  3. 13  shows  the  average  ion  energy  as  a  function  of  position  within  the 
discharge  at  the  beginning  (figure  3.13.a)  and  midway  through  (figure  3.13.b)  an  RF 
cycle.  Note  that  the  average  ion  energy  is  not  symmetric  during  the  cycle.  Figure  3. 13.a 
corresponds  to  a  time  just  after  the  sheath  at  the  left  electrode  has  fully  formed  and  the 
sheath  at  the  right  electrode  has  vanished.  In  this  manner,  the  ion  energy  responds  to  a 
half-cycle  average  -  during  the  half-cycle  when  a  sheath  is  formed  the  ions  in  that  sheath 
region  gain  more  energy  than  the  ions  at  the  opposite  end  of  the  discharge. 

Note  that  the  average  ion  energy  at  the  left  electrode  in  figure  3.13.a  is  approximately 
30  eV.  This  is  much  less  than  the  time  average  potential  available  to  the  ions  during  the 
preceding  half-cycle.  This  is  due  to  ion-neutral  collisions  which  preclude  most  of  the 
ions  from  reaching  the  wall  with  all  of  the  energy  available  from  the  sheath.  The  fast  ions 


in  the  sheath  collide  with  slow  neutrals  and  in  the  process  give  up  energy  to  the  neutrals. 
This  can  be  seen  in  figure  3. 14  which  is  a  z-vz  phase  space  plot  of  ions  in  the  sheath 
adjacent  to  the  left  electrode  at  the  beginning  of  an  RF  cycle.  The  ions  in  the  sheath  with 
the  most  energy  at  their  respective  positions  form  a  line  of  ions  which  have  not  yet 
undergone  a  collision.  If  one  extends  a  straight  line  through  this  line  of  ions  towards  the 
electrode,  the  line  intersects  the  electrode  at  about  9  x  104  to  105  m/sec.  This  corresponds 
to  roughly  175  to  210  V.  Thus,  the  ions  in  the  left  sheath  (at  least  at  the  beginning  of  the 
cycle)  must  have  responded  to  an  average  potential  of  around  200  V  (this  is  supported 
below  by  the  simulation  run  without  ion-neutral  collisions). 

The  average  number  of  collisions  an  ion  will  undergo  while  traversing  the  sheath,  C , 
may  be  calculated  given  the  collision  cross-section,  cr, ,  neutral  density,  Nn ,  and  the 

sheath  length,  d : 

C-dNnCfi  (3.9) 

Since  the  ion  collision  cross-section  used  in  this  simulation  is  independent  of  energy  (a 
constant  value  of  1.5  x  10'19  m2),  C  may  be  calculated  for  all  ions.  Given  the 
numerically  found  sheath  length  of  0.65  cm  (at  the  start  of  the  RF  cycle)  and  a  neutral 
density  of  8  x  1021  nr3  (0.25  torr)  this  gives  (on  average)  7.8  collisions  for  each  ion 
during  its  travel  across  the  sheath.  Thus,  the  average  ion  energy  at  the  wall  being  only 
15%  of  the  maximum  wall  energy  is  reasonable  given  this  degree  of  ion  collisionality  in 
the  sheath. 

Next,  the  results  from  a  simulation  with  the  same  conditions  as  above,  but  with  no 
ion-neutral  collisions,  are  presented.  Figure  3. 15  shows  the  average  ion  energy  at  the 
start  of  an  RF  cycle  as  a  function  of  position  within  the  discharge.  As  in  the  previous 
case,  the  energy  is  not  symmetric  about  center  of  the  discharge.  Unlike  the  previous  case, 
the  average  energy  at  the  left  electrode  is  222  eV  —  the  ions  don't  experience  any  energy 
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reducing  collisions  as  they  travel  in  the  sheath  towards  the  electrod  The  222  eV  ion 
energy  at  the  left  electrode  roughly  corresponds  to  the  maximum  value  of  ion  wall  energy 
for  the  case  with  ion  collisions.  The  z-vz  phase  space  plot  in  figure  3.16  corresponds  to 
the  beginning  of  an  RF  cycle  and  shows  the  ions  being  accelerated  through  the  sheath, 
experiencing  no  energy  reducing  collisions.  The  ions  which  have  not  attained  the 
maximum  energy  are  presumed  to  have  been  bom  in  the  sheath  as  a  result  of  electron- 
neutral  ionization. 

Since  the  ions  do  not  undergo  collisions  in  the  sheath  they  are  more  readily  lost  to  the 
electrodes.  This  results  in  a  lower  plasma  density  as  indicated  in  figure  3.17.  This 
decrease  in  plasma  density  has  implications  on  the  sheath  structure.  Figure  3. 18  shows 
the  electric  field  at  the  beginning  of  an  RF  cycle.  Note  that  the  sheath  length  at  the 
beginning  of  the  RF  cycle  has  increased  to  well  over  1  cm  as  compared  to  the  value  of 
0.65  cm  found  in  the  case  with  ion  collisions.  The  increased  sheath  length  with  the 
absence  of  ion  collisions  can  be  explained  as  follows.  Increased  ion  collisionality  in  the 
sheath  reduces  the  ion  energy  as  seen  in  figures  3.13  and  3.15.  In  order  to  maintain  the 
same  ion  flux,  the  ion  density  must  increase  for  increasing  ion  collisionality  (as  seen  in 
figure  3.1 1  and  3. 17).  This  increase  in  ion  density  for  increasing  ion  collisionality  results 


in  a  larger  electric  field  gradient,  — ~ ,  via  Poisson's  equation.  For  the  same  applied 

dx 


voltage  this  requires  a  smaller  sheath  length.  Therefore,  all  other  things  equal,  one 
expects  the  sheath  length  to  decrease  and  the  ion  density  to  increase  for  increasing  ion 
collisionality  (Sheridan  and  Goree,  1991:2801). 

Before  moving  on  to  the  two  dimensional  results,  the  results  from  a  parameter  study 
on  the  effect  of  varying  the  applied  voltage  is  presented.  The  experimental  conditions  are 
the  same  as  above  with  ion  collisions  included.  The  applied  voltage  is  varied  from  100  to 
1000  V.  The  plasma  density  at  the  center  of  the  discharge  and  the  sheath  length  are 
measured  as  a  function  of  applied  RF  voltage.  As  noted  above,  the  primary  factor  in 
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convergence  time  is  the  time  it  takes  to  achieve  the  correct  bulk  plasma  density  for  the 
given  discharge  conditions.  The  results  presented  here  can  be  used,  among  other  things, 
as  an  aid  in  choosing  the  initial  conditions.  Figure  3. 19  shows  the  plasma  density  in  the 
center  of  the  discharge  as  a  function  of  applied  RF  voltage  along  with  a  best  fit  line. 

Over  the  range  of  applied  voltages  the  plasma  density  is  approximately  proportional  to 
the  applied  voltage.  This  finding  does  not  agree  with  the  results  of  Vender  and  Boswell 
(see  equation  2.0)  for  their  simulation  of  atomic  hydrogen.  The  major  differences  in  the 
two  simulations  are  that  the  ionization  potential  of  hydrogen  is  approximately  10  eV  less 
than  that  of  helium,  and  they  did  not  include  elastic  collisions  in  their  model.  The  sheath 
length  at  the  beginning  of  an  RF  cycle  as  a  function  of  applied  voltage  is  shown  in  figure 
3.20.  The  sheath  length  is  fairly  constant  over  the  range  of  applied  voltages  for  all  but  the 
lowest  voltages  in  agreement  with  Vender  and  Boswell. 

Two  Dimensional  Results.  The  two-dimensional  model  was  run  with  the  same 
conditions  as  above  with  ion-neutral  charge  exchange  collisions.  Also,  the  discharge 
radius  is  set  to  2  cm.  Although  no  conclusions  can  be  drawn  from  this  data  due  to  the 
problems  sighted  above,  a  comparison  of  the  results  with  the  above  one -dimensional  data, 
as  well  as  investigating  the  ion  motion  in  the  sheaths,  is  worthwhile. 

Figure  3. 21. a  through  3.2 l.c  show  the  potential  as  a  function  of  position  in  the 
discharge  at  cat  /  2n  ~  0.25, 0.5,  and  0.75.  In  these  plots,  the  left  electrode  is  driven  and  is 
at  z=0  and  the  right  electrode  is  at  z=4  cm,  as  before.  The  discharge  tube  wall  is  at  r=2 
cm.  A  qualitative  comparison  to  the  previously  presented  one-dimensional  data  (see 
figures  3.9)  shows  a  similar  axial  profile  of  the  potential.  The  on-axis  noise  is  especially 
evident  at  tot  /  2rt  =  0.5.  The  radial  profile  of  the  potential  is  discussed  below. 

As  previously  discussed,  of  importance  to  RF  plasma-assisted  processing  is  the 
motion  of  the  ions  in  the  sheaths.  As  such,  the  remainder  of  this  section  is  devoted  to 
examining  the  ion  motion  in  the  sheaths  in  light  of  the  axial  and  radial  fields.  Figure  3.22 
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shows  the  average  ion  radial  velocity  as  a  function  of  position  within  the  discharge  at  the 
beginning  of  an  RF  cycle.  In  the  plasma  bulk  the  ions  have  small  radial  velocity 
components  (as  well  as  small  total  energies).  As  the  ions  enter  the  sheath  they  gain  radial 
velocity.  Upon  continuing  in  the  sheath  the  average  radial  velocity  decreases.  This  is 
presumed  to  be  due  to  the  increase  in  collision  frequency  with  increasing  velocity  -  as 
the  ions  gain  energy  from  the  sheath  they  are  more  likely  to  undergo  a  collision  thereby 
reducing  thei."  total  energy.  Since  the  amount  of  the  radial  potential  difference  in  the 
sheath  is  hard  to  determine  from  figures  3.21,  a  one-dimensional  plot  at  constant  z  =  3.5 
cm  at  tot  /  2ji  =  0.25  is  shown  in  figure  3.23.  This  corresponds  to  an  axial  position  just 
inside  the  sheath  at  this  time  in  the  RF  cycle.  One  can  see  a  large  potential  drop  in  going 
from  r  =  0  to  the  discharge  wall.  This  explains  the  increase  in  radial  velocity  at  the 
plasma-sheath  boundary.  Most  of  this  large  potential  drop  is  due  to  the  fact  that  the 
simulation  has  only  run  for  ten  RF  cycles  -  the  electrons  and  ions  are  still  diffusing  to  the 
glass  wall  at  independent  rates.  Unlike  the  electrons,  the  slower  ions  have  not  had  a 
chance  to  reach  the  glass  wall.  Therefore,  the  wall  is  more  negatively  charged  than  it 
would  be  at  some  later  time  after  the  ions  have  had  a  chance  to  set  up  a  space  charge 
electric  field  that  promotes  diffusion  of  both  species  at  the  am  bipolar  rate. 

The  average  ion  energy  as  a  function  of  axial  position  within  the  discharge  at  constant 
r  =  2.4  mm  at  the  beginning  of  an  RF  cycle  is  shown  in  figure  3.24.  Note  the  average  ion 
energy  at  the  left  electrode  of  jest  over  30  eV  and  the  asymmetry  similiar  to  the  one¬ 
dimensional  case  with  ion  collisions.  Figure  3.25  shows  the  average  ion  energy  across 
the  left  electrode.  Two-dimensional  data  like  this  is  required  to  determine  the  degree  of 
uniformity  of  ion  energy  deposition  across  the  wafer  platen. 

Figure  3.26  shows  the  average  ratio  of  the  radial  and  axial  velocity  components  for 
ions  at  the  left  electrode  (z  =  0)  at  the  beginning  of  an  RF  cycle.  This  plot  provides  data 
on  the  directionality  of  ions  at  the  electrode.  Note  that  for  values  of  radius  near  r  =  0  the 
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ions  have  velocities  directed  primarily  perpendicular  to  the  electrode.  As  the  glass  wall  is 
approached  through  increasing  values  of  r,  the  average  radial  velocity  component 
increases  to  approximately  18%  of  the  average  axial  component.  This  type  of  data  has 
important  implications  on  the  plasma  processing  characteristics  of  a  discharge. 

Further  evidence  of  the  anomalous  heating  due  to  the  noise  in  the  potential  is 
observed  through  the  average  electron  energy  in  the  center  of  the  discharge  in  the  two- 
dimensional  run.  The  observed  average  bulk  electron  energy  of  approximately  12  eV  is 
over  twice  as  large  as  the  value  of  5  eV  reported  in  the  above  one-dimensional  data  and 
by  other  researchers  under  similar  conditions. 

Summary  of  Results 

A  simple  test  case  involving  a  single  particle  was  compared  to  the  theoretical  solution 
with  good  results.  Furthermore,  interesting  observations  regarding  the  limitations  of 
particle  models  were  made.  The  single  particle  tests  give  confidence  that  major  portions 
of  the  model  (potential  solver,  field  solver,  particle  and  field  weighting,  and  portions  of 
the  particle  mover)  perform  correctly. 

A  pair  of  one-dimensional  discharge  runs  were  presented  that  provide  insight  into  the 
effect  of  ion  collisionality  on  the  characteristics  of  the  discharge's  sheath  structure.  A 
proportional  scaling  of  plasma  density  with  applied  RF  voltage  was  found  for  voltages 
over  the  range  of  100-1000  V  at  10  MHz  driving  frequency.  Problematic  two- 
dimensional  simulation  data  was  presented.  Although  no  conclusions  can  be  drawn  from 
this  data  it  provided  an  interesting  comparison  to  the  one-dimensional  data  as  well  as  an 
example  of  the  types  of  data  a  two-dimensional  simulation  of  RF  glow  discharges  makes 
available. 
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IV.  Suggestions  and  Recommendations 

This  chapter  is  a  list  of  required  and  suggested  improvements  to  the  model.  Each 
subsection  heading  is  phrased  as  a  suggestion  and  roughly  ordered  by  importance.  A 
basic  method  for  carrying  out  the  suggestion  is  given  in  the  text.  Also,  an  assessment  of 
the  level  of  difficulty  of  implementing  the  recommendation  is  provided. 

In  terms  of  general  suggestions  (not  specific  to  the  model  developed  for  this 
research),  a  fluid  model  should  be  obtained. 

Incorporate  a  Non-Uniform  Radial  Grid 

It  has  been  suggested  above  that  a  non-uniform  radial  grid  will  reduce  the  noise  in  the 
potential  on  the  radial  axis  of  the  grid.  By  forming  the  radial  grid  so  that  every  volume 
element  is  equal,  the  large  increase  in  volume  at  the  first  radial  grid  cell  will  be 
precluded. 

Since  the  current  subprogram  in  FISHPAK  will  only  handle  a  uniform  radial  grid, 
another  Poisson  solver  is  required.  It  has  been  found  that  a  Poisson  solver  that  uses  a 
non-uniform  radial  grid  can  be  constructed  by  using  some  of  the  lower  level  FISHPAK 
subprograms  and  a  coordinate  transformation.  The  incorporation  of  a  non-uniform  radial 
grid  also  requires  that  some  of  the  other  algorithms  be  revisited.  For  instance,  the  charge 
and  field  weighting  algorithms  will  require  modifications.  Also,  finite  differencing  the 
potential  to  find  the  radial  electric  field  components  becomes  a  bit  more  complicated 
since  each  grid  cell  does  not  have  the  same  radial  length. 

Before  incorporating  the  non-uniform  grid,  some  thought  should  be  given  to 
maintaining  the  relation  governing  the  choice  of  grid  spacing  in  equation  3.2.  If  this 
relation  is  held  for  the  first  grid  cell,  depending  upon  the  radius  of  the  simulated 
discharge  tube,  the  size  of  the  cells  near  the  glass  wall  become  increasingly  small  and  the 
total  number  of  cells  increases.  Furthermore,  the  accuracy  of  the  finite  differenced  radial 


electric  field  components  suffers  near  the  radial  axis.  Also,  the  accuracy  of  weighting  the 
charge  to  the  grid  and  the  field  components  to  the  charged  particles  suffers. 

As  a  quick  first  step  one  could  simulate  a  2D  cartesian  system.  Of  course,  curvature 
effects  are  lost  in  this  approach. 

Incorporate  Realistic  System  Boundaries 

The  current  model  incorporates  a  cylindrical  glass  wall  which  is  "butted  up  against" 
the  circular  electrodes.  This  geometry  was  chosen  primarily  based  on  the  simplicity  of 
ihe  resulting  boundary  conditions.  A  more  realistic  discharge  geometry  should  be 
adopted  along  with  the  option  for  a  grounded  conducting  wall  instead  of  the  currently 
implemented  glass  wall.  The  added  complexity  in  incorporating  the  more  realistic 
geometries  is  in  solving  Poisson's  equations  with  electrodes  that  are  not  on  the 
boundaries. 

Incorporate  Accurate  Modeling  of  Collision  Kinetics 

The  current  modeling  of  collision  kinetics  is  over  simplified.  Birdsall  outlines  a 
method  for  doing  the  energy  balance  and  finding  the  new  velocity  components  based  on 
an  approximation  to  the  differential  scattering  cross  section  (Birdsall,  1991:81).  This 
method  has  been  used  by  Surendra  and  Graves  and  is  suitable  for  noble  gases  (Surendra 
and  Graves,  1991:145)  Given  that  the  basic  structure  for  collision  processes  exists, 
inclusion  of  the  more  detailed  calculation  is  straightforward  and  should  be  incorporated. 

Optimize  the  Code 

Work  should  be  done  to  optimize  the  code  for  run-time.  As  the  model  was  coded, 
run-time  was  a  secondary  consideration.  Emphasis  was  put  on  ensuring  the  code  was 
correct  and  readable.  So  that  longer  computer  runs  can  be  made,  the  code  should  execute 
as  fast  as  possible.  Fust,  a  code  profiler  should  be  used  to  determine  the  portions  of  the 
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code  were  a  majority  of  computation  time  is  spent.  These  portions  of  the  code  should 
then  be  considered  for  run-time  optimization. 

It  was  noted  above  that  the  C  to  FORTRAN  interface  requires  that  the  two- 
dimensional  array  used  to  store  the  grid  quantities  be  converted  from  row  major  format  to 
column  major  format  and  back  to  row  major  format  each  time-step.  This  is  an  immediate 
candidate  for  optimization.  Unfortunately,  the  remedy  for  this  situation  requires  that 
either  the  model  be  recoded  in  FORTRAN  or  the  Poisson  solver  be  recoded  using  C. 
Before  either  of  these  alternatives  are  done,  the  relative  amount  of  time  spent  converting 
between  storage  formats  should  be  determined.  It  may  tum  out  that  it  is  a  small  fraction 
of  the  total  time  in  the  time-step  loop  and  optimizing  other  portions  of  the  code  may  have 
more  pay-off. 

Port  the  Model  to  a  Faster  Computer  with  Large  Main  Memory 

Since  the  software  does  not  take  advantage  of  any  system  specific  features,  porting 
the  code  to  other  computer  systems  will  be  straightforward.  Depending  upon  the 
computing  resource  available,  this  may  allow  simulation  runs  involving  more  particles 
and  a  more  accurate  grid  spacing.  Furthermore,  once  the  above  mentioned  problems  are 
corrected,  quantitative  comparisons  to  previous  work  and  theory  can  be  made  for  fully 
converged  2D  runs. 

Include  the  Effects  of  an  Externally  Applied  Magnetic  Field 

External  magnetic  fields  are  used  in  some  devices  to  confine  the  plasma.  These 
effects  can  be  included  by  modifying  the  particle  mover  to  include  the  magnetic  field 
term  in  the  Newton-Lorentz  equation  of  motion. 

Add  External  Circuit  Elements 

Currently,  the  simulation  does  not  include  a  model  of  an  external  circuit.  The 
complexity  of  adding  an  external  circuit  simulation  to  the  particle  model  has  not  been 
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investigated  and  its  impact  on  computer  run-time  is  not  yet  known.  Given  this 
uncertainty,  more  research  into  this  area  is  required  before  an  external  circuit  is  added  to 
the  existing  model.  As  a  first  step,  Lawson  provides  a  general  outline  of  adding  an 
external  RLC  circuit  to  a  particle  simulation  (Lawson,  1989:267-272). 
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Appendix  A:  Electric  Field  Within  a  Semi-Infinite  Dielectric 


In  this  appendix  it  will  be  shown  that  the  normal  component  of  the  electric  field 
within  a  semi-infinite  dielectric  (due  to  a  point  charge  opposite  the  dielectric)  becomes 
vanishingly  small  for  large  dielectric  constant.  Using  the  method  of  images  and 
following  the  development  of  Becker  and  Sommerfeld,  consider  the  configuration  shown 
in  figure  A.l  (Becker,  1964:99-100;  Sommerfeld,  1964:63-64).  Shown  is  a  particle  of 
charge  e  in  air  (region  1  of  dielectric  constant  e,)  opposite  a  semi-infinite  dielectric 
(region  2  of  dielectric  constant  e2 )  along  with  images  of  charge  -e'  and  e" .  The 
potential  within  the  dielectric  is  obtained  by  considering  only  image  charge  e"  at  point  A 
and  as  if  the  dielectric  occupied  all  space  (Becker,  1964:99): 


4  Jt  e2  r 


(A.l) 


The  potential  in  region  1  is  due  to  the  real  charge  at  A  and  the  image  charge  at  B  (Becker, 
1964:99): 

e  e' 

- - - r  (A-2) 

4  jt  e,  r  4  ;r  e,  r 

By  applying  the  appropriate  boundary  conditions  expressions  for  the  image  charges 
may  be  obtained.  First,  the  boundary  condition  on  the  normal  component  of  the  electric 
displacement  vector,  D.  Let  x  -  0  be  at  the  interface  and  Jf  /  /  -  h  so  that  the  boundary 
condition  on  D  at  x  -0  becomes  (Sommerfeld,  1964:63): 


1  dx  !  ax 

Carry  out  the  following  derivatives  (Sommerfeld,  1964:63): 


(A  .3) 


x~a 
r3 


& 


dx\r'J 


x  +a 


(A.4) 


so  that  equation  A3  can  be  rewritten  as: 
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{x-  a)  e  (j x+a)e'  ( x-a)e " 


ft 


(A.5) 


r  r  r 

If  equation  A.5  is  evaluated  at  x  «  0  one  obtains  the  following  expression  (Becker, 
1964:100): 


e  +  e  - e 

Next,  apply  the  boundary  condition  on  the  potential,  namely 

-  <f>2  at  x  -  0 

which  gives  the  following  expression  (Becker,  1964: 100): 

e-e'„h.e" 


(A. 6) 


(A- 7) 


(A  .8) 


Equations  A.6  and  A.8  can  then  be  used  to  find  the  image  charge  at  point  A  in  terms  of 
the  real  charge  at  point  A  and  the  dielectric  constants  (Becker,  1964: 100): 

2  £, 


e  -  e- 


e,+e2 


(A-9) 


Equation  A.9  can  now  be  used  to  find  the  normal  component  of  the  electric  field  within 
the  dielectric.  For  x  <  0 : 

(<»-*) 


(«i+«a)' 


(A.  10) 


In  the  limit  of  region  2  having  infinite  dielectric  constant  the  normal  component  of  the 
electric  field  within  region  2  vanishes. 
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Appendix  B:  Sample  Input  Data  File 


dz 

2e-4 

Nz 

200 

dr 

2e-4 

Nr 

75 

dt 

le-10 

Nt 

10000 

npr 

6 

npz 

6 

MiMe 

7344.0 

Ve 

1.0e6 

Vi 

2 . 5e3 

P 

300e-3 

Ns 

8 . 0e4 

sigma_total 

7.5e-14 

sigma_max_ion 

3 . 5e-21 

eOion 

24.6 

elion 

70.0 

e2ion 

110.0 

sigmamaxelas 

5.5e-20 

eOelas 

4.0 

sigma_ch_ex 

1.5e-19 

SecE 

0.0 

Sec  I 

0.3 

VSec 

5.93e5 

RefE 

0.25 

RefI 

0.0 

VrfO 

500.0 

frf  0 

le7 

VdcO 

0.0 

VrfL 

0.0 

frfL 

0.0 

VdcL 

0.0 

weighting 

2 

Nout 

50 

Pout 

10000 

Eout 

10000 

lout 

10000 

Fout 

10000 

Sout 

10000 

Cout 

0 

Rout 

0 

Tout 

0 

seed 

12556789 
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Appendix  C:  Diagrams 


1 

y  Pumping 
System 


Figure  1.1 .  Schematic  of  a  typical  RF  glow  discharge  (Flamm  and  Herb,  1989: 15). 
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Mask  Etch 


Substrate 


Figure  1.2m.  Characteristics  of  the  isotropic,  wet  etched  pattern  (National  Research 
Council,  1991:36). 


Figure  1.2.b.  Characteristics  of  the  anisotropic,  dry  etched  pattern  (National  Research 
Council,  1991:36). 
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AND  WAFER  PLATEN 


Figure  U.  Schematic  of  the  "Reinberg  reactor."  Note  the  radial  gas  flow  feature.  This 
design  was  patented  by  Reinberg  in  1975  (Hamm  and  Herb,  1989:66). 
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Neutral  Free  Radical  Etching 


Ion-Assisted  Etching 


Neutral 


Figure  1.4. 


Sputter  Etching 


Basic  plasma  assisted  etching  mechanisms  (Flamm  and  Herb,  1989: 19). 
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Integrate  equations  of  motion 


Figure  2.1.  Computational  cycle  used  in  particle  simulations  with  Monte  Carlo 
collisions  (Birdsall,  1990:81). 
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Figure  2.3.  Volume  elements  used  in  forming  the  charge  density.  The  charge  density  is 
weighted  to  the  dots  centered  on  the  column  elements  at  the  grid  points.  The  volume 
elements  associated  with  on-axis  grid  points  have  a  different  form  than  the  off-axis 
elements  (Birdsall  and  Langsdon,  1990:333). 
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Figure  2.4.  Pictorial  representation  of  the  leapfrog  integration  scheme.  Note  that  the 
forces  and  positions  are  known  at  integral  time-steps  while  the  velocities  are  known  at 
half-integral  time-steps  (Birdsall  and  Langsdon,  1991:13). 
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Figure  2.5.  Grid  with  boundaries  shown.  The  boundary  at  z  =  0  corresponds  to  the  left 
electrode  while  the  boundary  at  z  =  L  corresponds  to  the  right  electrode.  The  boundary  at 
r  =  R  is  the  glass  wall  of  the  discharge  tube. 


Surface 

Charge 


Figure  2.7.  Diagram  used  to  explain  the  bilinear  and  volume  weighting  schemes 
(Biidsall  and  Langdon,  1991:309). 
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Figure  2.9.  Form  of  elastic  collision  cross-section  with  parameters  discussed  in  text 
(Birdsall,  1991:81). 
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Figure  3.1.  Location  of  electron,  grounded  electrodes  and  first  four  image  charges  for 
single  particle  test  case  1. 
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Figure  3.2.  Location  of  an  electron,  grounded  electrodes,  and  six  image  charges  for 
single  particle  test  case  2.  The  distance  beneath  each  image  charge  indicatestance  of  the 
image  from  the  electron  at  z. 
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Figure  33.  Percent  difference  between  analytically  calculated  and  model  calculated 
acceleration  vs  simulation  time-step  for  single  particle  test  case  2  using  bilinear 
weighting. 


Normalised  z  (unitless) 


Figure  3.4.  Normalized  z-component  of  position  vs  simulation  time-step  for  single 
particle  test  case  2  using  bilinear  weighting. 
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Figure  3.5.  Percent  difference  between  analytically  calculated  and  model  calculated 
acceleration  vs  simulation  time-step  for  single  particle  test  case  2  using  NGP  weighting. 
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Figure  3.6.  Normalized  z-component  of  position  vs  simulation  time-step  for  single 
particle  test  case  2  using  NGP  weighting. 
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Figure  3.7 .a.  Potential  at  the  start  of  the  eleventh  RF  cycle.  Note  the  on-axis  noise. 


Figure  3.7.b.  One  dimensional  slice  at  constant  radius  (r  =  0)  of  potential  shown  in 
figure  3.7.a. 
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1000  Particles  Per  Grid  Cell 


Figure  3.8.c.  Potential  for  conditions  discussed  in  chapter  three  with  1000  particles  per 
grid  cell. 


1000  Particles  Per  Grid  Cell 


Figure  3.8.d.  Potential  for  particles  placed  off-axis  with  1000  particles  per  grid  cell . 


cot  /  2n  =  0.250 


Figure  3.9 .a.  Potential  as  a  function  of  position  within  the  discharge  at  three  times 
between  the  beginning  and  one-quarter  of  the  way  through  the  RF  cycle. 


oit  /  2ji  =  0.50 


Figure  3.9.b.  Potential  as  a  function  of  position  within  the  discharge  at  three  times 
between  one-half  and  three-quarters  of  the  way  through  the  RF  cycle. 
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Figure  3.10.a.  Electric  field  as  a  function  of  position  within  the  discharge  at  three  times 
between  the  beginning  and  one-quarter  of  the  way  through  the  RF  cycle . 


Figure  3.10.b.  Electric  field  as  a  function  of  position  within  the  discharge  at  three  times 
between  one-half  and  three-quarters  of  the  way  through  the  RF  cycle . 
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Figure  3.1  l.a.  Charged  particle  densities  as  a  function  of  position  in  the  discharge  at  the 
beginning  of  an  RF  cycle. 


Figure  3.1  l.b.  Charged  particle  densities  as  a  function  of  position  in  the  discharge  one- 
eighth  of  the  way  through  an  RF  cycle. 
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Figure  3.1  l.c.  Charged  particle  densities  as  a  function  of  position  in  the  discharge  one- 
quarter  of  the  way  through  an  RF  cycle. 
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Figure  3.1241.  Average  electron  energy  as  a  function  of  position  within  the  discharge  at 
cot  /  2ji  =  0.0. 


Figure  3.12.b.  Average  electron  energy  as  a  function  of  position  within  the  discharge  at 
cot /2ji  =  0.125. 
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Figure  3.12.C.  Average  electron  energy  as  a  function  of  position  within  the  discharge  at 
cot  /  2rt  =  0.25. 


Figure  3.12.d.  Average  electron  energy  as  a  function  of  position  within  the  discharge  at 
cut  /  2ji  =  0.375. 
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Figure  3.12.e.  Average  electron  energy  as  a  function  of  position  within  the  discharge  at 
cat  /  2ji  =  0.5. 
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Figure  3.13.a.  Average  ion  energy  as  a  function  of  position  within  the  discharge  at 
beginning  of  RF  cycle. 


Figure  3.13.b.  Average  ion  energy  as  a  function  of  position  within  the  discharge  midway 
through  the  RF  cycle. 
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Figure  3.15.  Average  ion  energy  at  the  beginning  of  the  RF  cycle  as  a  function  of 
position  within  the  discharge  for  the  case  with  no  ion  collisions. 
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Figure  3.16.  Ion  z-vz  phase  space  at  the  beginning  of  the  RF  cycle  for  the  case  with  no 
ion  collisions  for  the  left  half  of  the  discharge  (Om  s  zi  0.02  m). 
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Figure  3.18.  Electric  field  at  the  beginning  of  the  RF  cycle  for  the  case  with  no  ion 
collisions. 
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Figure  3.21.a.  Potential  as  a 
for  the  two-dimensional  case. 


Figure  3.21.b.  Potential  as  a 
the  two-dimensional  case. 


Figure  3.21.C.  Potential  as  a  function  of  position  within  the  discharge  at  cot  /  2jc  =  0.75 
for  the  two-dimensional  case. 
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Figure  3.22.  Average  ion  radial  velocity  as  a  function  of  position  within  the  discharge  at 
the  beginning  of  an  RF  cycle  for  the  two-dimensional  case. 
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Average  Ion  Energy  (eV) 


Figure  3.24.  Average  ion  energy  as  a  function  of  axial  position  within  the  discharge  at 
constant  r  =  2.4  mm  at  the  beginning  of  an  RF  cycle  for  the  two-dimensional  case. 
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Figure  3-25.  Average  ion  energy  across  the  left  electrode  for  the  two-dimensional  case. 
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Figure  3.26.  Average  ratio  of  the  radial  and  axial  velocity  components  for  ions  at  the  left 
electrode  (z  =  0)  at  the  beginning  of  an  RF  cycle  for  the  two-dimensional  case. 
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Figure  A.l.  Diagram  referred  to  in  appendix  A  (Becker,  1964:99-100). 
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